function [varargout] = cub_square_padua(n,xyrange,varargin)

%--------------------------------------------------------------------------
% USAGE of "pdcub".
%
% PadL = pdcub(n,xyrange)
% [cubature,PadL] = pdcub(n,xyrange,funct,opt1,opt2,...)
%
% This function computes the Padua points defined in the rectangle 
% [xyrange(1),xyrange(2)] x [xyrange(3),xyrange(4)] and the 
% corresponding cubature weights as a matrix of abscissas (first
% column), ordinates (second column) and weights (third column).
% Optionally, it computes the cubature of the function f (optional 
% input argument).
%--------------------------------------------------------------------------
% INPUT.
%
% n        : interpolation degree
% xyrange  : a vector [a,b,c,d] defining the rectangle [a,b] x [c,d]
% funct    : function to be interpolated in the form 
%            funct(x,y,opt1,opt2,...), where opt1, opt2, ... are
%            optional arguments for funct 
%
% OUTPUT.
%
% PadL     : a matrix with the abscissas of Padua points in the 
%            first column, the ordinates in the second and the 
%            cubature weights in the third
% cubature : cubature of the integrand funct
%--------------------------------------------------------------------------
% EXAMPLES.
%
% 1)
% Compute the Padua points and the cubature weights of degree 50
%
% xyrange = [0,1,0,1];
% PadL = pdcub(50,xyrange);
%
% 2)
% Compute the cubature of the Franke's function defined in funct.m
%
% xyrange = [0,1,0,1];
% [cubature,PadL] = pdcub(50,xyrange,@f);
%--------------------------------------------------------------------------
% FUNCTIONS CALLED BY THIS CODE:
% 1. pdpts
% 2. pdwtsMM
%--------------------------------------------------------------------------
% Copyright (C) 2008-2009 
% Marco Caliari, Stefano De Marchi, 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 2 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:  
%          Marco Caliari     <marco.caliari@univr.it>
%          Stefano De Marchi <stefano.demarchi@univr.i>   
%          Alvise Sommariva  <alvise@euler.math.unipd.it>
%          Marco Vianello    <marcov@euler.math.unipd.it>   
%
% Date: March 14, 2009.
% Modified: November 4, 2025 (A. Sommariva)
%--------------------------------------------------------------------------

% Compute the Padua points in the rectangle defined by xyrange
PadL = pdpts(n,xyrange);
% Compute the cubature weights
PadL(:,3) = pdwtsMM(n,xyrange);


if (nargin < 2)
  error('Too few input arguments')
elseif (nargin >= 3)
  funct = varargin{1};
  varargout{1} = PadL(:,3)'*feval(funct,PadL(:,1),PadL(:,2),...
      varargin{2:end});
  varargout{2} = PadL;
else
  varargout{1} = PadL;
end

%--------------------------------------------------------------------------
% OCTAVE TESTS.
%--------------------------------------------------------------------------
% Octave testing: type
%
% test pdcub
%
% at the Octave prompt
%
%!test
%! disp('Degree 0 cubature')
%! xyrange = [0,1,0,1];
%! [cubature,PadL] = pdcub(0,xyrange,@funct,13);
%! expected = 1;
%! assert(cubature/expected,1,10*eps);
%!test
%! disp('Franke''s function cubature')
%! xyrange = [0,1,0,1];
%! [cubature,PadL] = pdcub(500,xyrange,@funct);
%! expected = 4.069695894915573e-01;
%! assert(cubature/expected,1,20*eps);  




function [varargout] = pdpts(n,varargin)

%--------------------------------------------------------------------------
% USAGE of "pdpts".
%
% Pad = pdpts(n)
% Pad = pdpts(n,xyrange)
% [X1,Y1,X2,Y2] = pdpts(n)
% [X1,Y1,X2,Y2] = pdpts(n,xyrange)
%
% Compute the (first family of) Padua points, either as a matrix Pad 
% with their abscissas in the first column and their ordinates in the
% second, or as two subgrids X1,Y1 and X2,Y2, respectively.
%--------------------------------------------------------------------------
% INPUT.    
%
% n           : interpolation degree
% xyrange     : an optional vector [a,b,c,d] defining the rectangle 
%               [a,b] x [c,d]. Otherwise, xyrange = [-1,1,-1,1]
%
% OUTPUT.  
%
% Pad         : matrix of size ((n+1)*(n+2)/2) x 2 such
%               that (Pad(:,1),Pad(:,2)) defines the Padua points in the
%               rectangle [xyrange(1),xyrange(2)] x [xyrange(3),xyrange(4)].  
% X1,Y1,X2,Y2 : the two subgrids X1,Y1 and X2,Y2 defining the Padua points
%--------------------------------------------------------------------------
% FUNCTIONS CALLED BY THIS CODE:
% no external function is used.
%--------------------------------------------------------------------------
% Copyright (C) 2008-2009 
% Marco Caliari, Stefano De Marchi, 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 2 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:  
%          Marco Caliari     <marco.caliari@univr.it>
%          Stefano De Marchi <stefano.demarchi@univr.i>   
%          Alvise Sommariva  <alvise@euler.math.unipd.it>
%          Marco Vianello    <marcov@euler.math.unipd.it>   
%
% Date: March 09, 2009.
%--------------------------------------------------------------------------

if (nargin == 1)
% standard square [-1,1] x [-1,1]
  xyrange = [-1,1,-1,1];
else
% rectangle [xyrange(1),xyrange(2)] x [xyrange(3),xyrange(4)]
  xyrange = varargin{1};
end
if (n == 0)
% degree 0
  if (nargout ~= 4)
% points as a single matrix
    Pad = [xyrange(1),xyrange(3)];
    varargout = {Pad};
  else
% points as two subgrids
    X1 = xyrange(1);
    Y1 = xyrange(3);
    X2 = zeros(1,0);
    Y2 = zeros(1,0);
    varargout = {X1,Y1,X2,Y2};
  end  
else
% degree > 0
  zn = (xyrange(1)+xyrange(2)+(xyrange(2)-xyrange(1))*...
       cos(linspace(0,1,n+1)*pi))/2;
  zn1 = (xyrange(3)+xyrange(4)+(xyrange(4)-xyrange(3))*...
	cos(linspace(0,1,n+2)*pi))/2;
  if (nargout ~= 4)
% points as a single matrix
    [Pad1,Pad2] = meshgrid(zn,zn1);
    [M1,M2] = meshgrid([0:n],[0:n+1]);
    findM = find(mod(M1+M2,2));
    Pad = [Pad1(findM),Pad2(findM)];
    varargout = {Pad};
  else
% points as two (mesh)grids  
    En = zn(1:2:n+1);
    On = zn(2:2:n+1);
    En1 = zn1(1:2:n+2);
    On1 = zn1(2:2:n+2);
    [X1,Y1] = meshgrid(En,On1);
    [X2,Y2] = meshgrid(On,En1);
    varargout = {X1,Y1,X2,Y2};
  end  
end

%--------------------------------------------------------------------------
% OCTAVE TESTS.
%--------------------------------------------------------------------------
% Octave testing: type
%
% test pdpts
%
% at the Octave prompt
%
%!test
%! disp('Degree 0')
%! Pad = pdpts(0);
%! expected = [-1,-1];
%! assert(Pad,expected,10*eps); 
%!test
%! disp('Degree 1')
%! Pad = pdpts(1);
%! expected = [cos([0;1;1]*pi),cos([1;0;2]*pi/2)];
%! assert(Pad,expected,10*eps); 
%!test
%! disp('Degree 2 with xyrange')
%! Pad = pdpts(2,[0,1,0,2]);
%! expected = [(cos([0;0;1;1;2;2]*pi/2)+1)/2,cos([1;3;0;2;1;3]*pi/3)+1];
%! assert(Pad,expected,10*eps); 






function [varargout] = pdwtsMM(n,varargin)

%--------------------------------------------------------------------------
% USAGE of "pdwtsMM".
%
% L = pdwtsMM(n)
% L = pdwtsMM(n,xyrange)
% [L1,L2,L] = pdwtsMM(n)
% [L1,L2,L] = pdwtsMM(n,xyrange)
%
% Compute the cubature weights L so that, if Pad is the
% matrix of Padua points computed through the call
% Pad = pdpts(n) or Pad = pdpts(n,xyrange), then the cubature of 
% the function f is given by L'*f(Pad(:,1),Pad(:,2)). 
% Otherwise, one can compute the cubature weights L1 and L2, so 
% that, if X1,Y1 and X2,Y2 are the subgrids of Padua points 
% computed through the call [X1,Y1,X2,Y2] = pdpts(n) or 
% [X1,Y1,X2,Y2] = pdpts(n,xyrange), then the cubature of the 
% function f is given by 
% sum(sum(L1.*f(X1,Y1)))+sum(sum(L1.*f(X2,Y2))).
%--------------------------------------------------------------------------
% INPUT.    
%
% n       : interpolation degree
% xyrange : an optional vector [a,b,c,d] defining the rectangle 
%           [a,b] x [c,d]. By default, xyrange = [-1,1,-1,1]
%
% OUTPUT.   
%
% L       : cubature weights associated to the matrix Pad of 
%           Padua points
% L1,L2   : cubature weights associated to the subgrids 
%           X1,Y1 and X2,Y2 of Padua poins.
%--------------------------------------------------------------------------
% FUNCTIONS CALLED BY THIS CODE:
% no external function is used.
%--------------------------------------------------------------------------
% Copyright (C) 2008-2009 
% Marco Caliari, Stefano De Marchi, 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 2 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:  
%          Marco Caliari     <marco.caliari@univr.it>
%          Stefano De Marchi <stefano.demarchi@univr.i>   
%          Alvise Sommariva  <alvise@euler.math.unipd.it>
%          Marco Vianello    <marcov@euler.math.unipd.it>   
%
% Date: March 10, 2009.
%--------------------------------------------------------------------------

if (nargin == 1)
  xyrange = [-1,1,-1,1];
else
  xyrange = varargin{1};
end
if (n == 0)
% degree 0
  L1 = (xyrange(2)-xyrange(1))*(xyrange(4)-xyrange(3));
  L2 = zeros(1,0);
  L = L1;
else
  argn = linspace(0,pi,n+1);
  argn1 = linspace(0,pi,n+2);
  k = [0:2:n]';
  l = (n-mod(n,2))/2+1;
% even-degree Chebyshev polynomials on the subgrids
  TE1 = cos(k*argn(1:2:n+1));
  TE1(2:l,:) = TE1(2:l,:)*sqrt(2);
  TO1 = cos(k*argn(2:2:n+1));
  TO1(2:l,:) = TO1(2:l,:)*sqrt(2);
  TE2 = cos(k*argn1(1:2:n+2));
  TE2(2:l,:) = TE2(2:l,:)*sqrt(2);
  TO2 = cos(k*argn1(2:2:n+2));
  TO2(2:l,:) = TO2(2:l,:)*sqrt(2);
% even,even moments matrix
  mom = 2*sqrt(2)./(1-k.^2);
  mom(1) = 2;
  [M1,M2] = meshgrid(mom);
  M = M1.*M2;
  Mmom = fliplr(triu(fliplr(M)));
% interpolation weights matrices
  W1 = 2*ones(l)/(n*(n+1));
  W2 = 2*ones((n+mod(n,2))/2+1,(n+mod(n,2))/2)/(n*(n+1));
  W1(:,1) = W1(:,1)/2;
  W2(1,:) = W2(1,:)/2;
  if (mod(n,2) == 0)
    Mmom(n/2+1,1) = Mmom(n/2+1,1)/2;
    W1(:,n/2+1) = W1(:,n/2+1)/2;
    W1(n/2+1,:) = W1(n/2+1,:)/2;
  else
    W2((n+1)/2+1,:) = W2((n+1)/2+1,:)/2;
    W2(:,(n+1)/2) = W2(:,(n+1)/2)/2;
  end
% cubature weights as matrices on the subgrids.
  L1 = W1.*(TE1'*Mmom*TO2)';
  L2 = W2.*(TO1'*Mmom*TE2)';
  if (mod(n,2) == 0)
    L = zeros(n/2+1,n+1);
    L(:,1:2:n+1) = L1;
    L(:,2:2:n+1) = L2;
    L = L(:);
  else
    L = zeros((n+1)/2,(n+2));
    L = [L1',L2']';
    L = L(:);
  end
  L = L*(xyrange(2)-xyrange(1))*(xyrange(4)-xyrange(3))/4;
end
if (nargout == 0 | nargout == 1)
  varargout{1} = L;
else
  varargout{1} = L1;
  varargout{2} = L2;
  varargout{3} = L;
end

%-------------------------------------------------------------------------------
% OCTAVE TESTS.
%-------------------------------------------------------------------------------
% Octave testing: type
%
% test pdwtsMM
%
% at the Octave prompt
%
%!test
%! disp('Degree 0 weight')
%! xyrange = [0,1,0,1];
%! L = pdwtsMM(0,xyrange);
%! expected = 1;
%! assert(L/expected,1,10*eps);
%!test
%! disp('Degree 0 weight, subgrids')
%! xyrange = [-1,1,-1,1];
%! [L{1},L{2}] = pdwtsMM(0,xyrange);
%! expected{1} = 4; 
%! expected{2} = zeros(1,0);
%! assert(L,expected,10*eps);
%!test
%! disp('Degree 1 weights')
%! xyrange = [0,1,0,1];
%! L = pdwtsMM(1,xyrange);
%! expected = [0.5;0.25;0.25];
%! assert(L,expected,10*eps);
%!test
%! disp('Degree 1 weights, subgrids')
%! xyrange = [-1,1,-1,1];
%! [L{1},L{2}] = pdwtsMM(1,xyrange);
%! expected{1} = 2;
%! expected{2} = [1;1];
%! assert(L,expected,10*eps);
%!test
%! disp('Degree 2 weights')
%! xyrange = [0,1,0,1];
%! L = pdwtsMM(2,xyrange);
%! expected = [1/6;0;1/9;5/9;1/6;0];
%! assert(L,expected,10*eps);
%!test
%! disp('Degree 2 weights, subgrids')
%! xyrange = [-1,1,-1,1];
%! [L{1},L{2}] = pdwtsMM(2,xyrange);
%! expected{1} = [2/3,2/3;0,0];
%! expected{2} = [4/9;20/9];
%! assert(L,expected,10*eps);
