
function demo_cubature_2009vs2019(n)

% OBJECT:
% the procedure tests the cubature on curvilinear polygons

% ......................... initial data ..................................

if nargin < 1
n=6; % degree of precision
end
extraction_type=0; % 1: lsqnonneg, 2: fast LH, 0: backslash
example=14;
function_type=4;
plot=1;

SPLtypestring='not-a-knot';

% ......................... comparison ....................................
polygon_sides=define_polygon(example);
spline_order_vett=[4 size(polygon_sides,1)];

astr=num2str(rand(1),15); bstr=num2str(rand(1),15);
nstr=num2str(n,15);
fstr=strcat('(',astr,'*x+',bstr,'*y).^',nstr);
f=inline(fstr);

% ..... splinegauss_2009b .....
tic;
[nodes_x,nodes_y,weights]=splinegauss_2009b(n,polygon_sides,...
    [],[],[],spline_order_vett,0,SPLtypestring,4);
tt(1)=toc;
fnodes=feval(f,nodes_x,nodes_y);
I_2009=weights'*fnodes;

% ..... splcub .....
X=polygon_sides(:,1); Y=polygon_sides(:,2);
tic;
[Sx,Sy]=Sx_Sy_preparation(X,Y,1,spline_order_vett,SPLtypestring);
[xyw,res,Z,Zin] = splcub(n,Sx,Sy,extraction_type);
tt(2)=toc;
fnodes=feval(f,xyw(:,1),xyw(:,2));
I_2019=xyw(:,3)'*fnodes;

% ..... compare results .....
fprintf('\n \t   f(x)=')
fprintf(fstr);
fprintf('\n \t 2009: %1.15e',I_2009);
fprintf('\n \t 2019: %1.15e',I_2019);
fprintf('\n \t rel.err.: %1.15e',abs(I_2009-I_2019)/abs(I_2019));

fprintf('\n \n \t ...... properties ......');
fprintf('\n \t ** ade: %3.0f',n);
fprintf('\n \t ** compression residual: %1.15e',res);
cond2009=sum(abs(weights))/sum(weights);
cond2019=sum(abs(xyw(:,3)))/sum(xyw(:,3));
fprintf('\n \t * cub.cond. 2009: %1.3e',cond2009);
fprintf('\n \t   cub.cond. 2019: %1.3e',cond2019);
fprintf('\n \t * cputime 2009: %1.1e',tt(1));
fprintf('\n \t   cputime 2019: %1.1e',tt(2));
fprintf('\n \t   ratio cputime: %1.1e',tt(2)/tt(1));
fprintf('\n \t * cardinality 2009: %10.0f',length(nodes_x));
fprintf('\n \t   cardinality 2019: %10.0f',size(xyw,1));
in_2009=incurvpolygon([nodes_x nodes_y],Sx,Sy);
in_2019=incurvpolygon(xyw(:,1:2),Sx,Sy);
fprintf('\n \t * internal pts 2009: %9.0f',length(find(in_2009 == 1)));
fprintf('\n \t   internal pts 2019: %9.0f',length(find(in_2019 == 1)));
fprintf('\n \t * spl 2019 all pts: %9.0f',size(Z,1));
fprintf('\n \t   spl 2019 in  pts: %9.0f',size(Zin,1));
fprintf('\n \n');

% ..... plot domain and pointsets .....
plot_spl(Sx,Sy,xyw(:,1:2),Z,Zin);






function [Sx,Sy,Nsub]=Sx_Sy_preparation(X,Y,Nsub,spline_parms,SPLtypestring)

% OBJECT:
%
%
% INPUT:
% X,Y:
% Nsub:
% spline_parms:
% SPLtypestring:
%
% OUTPUT:
% Sx,Sy:
% Nsub:
%

% check input variables.
L=size(spline_parms,1);

for i=1:L % define spline
    
    if i == 1
        imin=1;
    else
        imin=spline_parms(i-1,2);
    end
    imax=spline_parms(i,2);
    
    tL=imin:imax;
    xL=X(imin:imax);
    yL=Y(imin:imax);
    
    [SxL,SyL]=compute_parametric_spline(tL,xL,yL,...
        spline_parms(i,1),SPLtypestring);
    
    Sx(i)=SxL;
    Sy(i)=SyL;
    
end





%--------------------------------------------------------------------------
% compute_parametric_spline
%--------------------------------------------------------------------------

function [ppx,ppy]=compute_parametric_spline(s,x,y,spline_order,...
    SPLtypestring)

% OBJECT:
% compute parametric spline relavant parameters "ppx", "ppy" so that a
% point at the boundary of the  domain has coordinates (x(s),y(s))
%
% INPUT:
% s: parameter data.
% x: determine spline x(s) interpolating (s,x)
% y: determine spline y(s) interpolating (s,y)
% spline_order: spline order (i.e. degree + 1)
% SPLtypestring: string with the spline type i.e.
%             'complete'   : match endslopes (as given in VALCONDS, with
%                     default as under *default*).
%             'not-a-knot' : make spline C^3 across first and last interior
%                     break (ignoring VALCONDS if given).
%             'periodic'   : match first and second derivatives at first
%                     data point with those at last data point (ignoring
%                     VALCONDS if given).
%             'second'     : match end second derivatives (as given in
%                    VALCONDS, with default [0 0], i.e., as in variational).
%             'variational': set end second derivatives equal to zero
%                     (ignoring VALCONDS if given).
%
% OUTPUT:
% ppx: spline x(t) data
% ppy: spline y(t) data


switch spline_order
    case 4
        
        % CUBIC SPLINES BY CSAPE. AS DEFAULT WE USE PERIODIC CUBIC SPLINES.
        % Derivatives parameters are computed as well.
        ppx=csape(s,x,SPLtypestring);
        ppy=csape(s,y,SPLtypestring);
        
    otherwise
        
        ppx=spapi(spline_order,s,x);
        ppy=spapi(spline_order,s,y);
        
end

if (ppx.form =='B-')
    ppx=fn2fm(ppx,'pp');
    ppy=fn2fm(ppy,'pp');
end





function polygon_sides=define_polygon(polygon_type)

%-------------------------------------------------------------------------------
% POLYGON DEFINITION.
%-------------------------------------------------------------------------------
%
% INPUT:
%
% [polygon_type]: PARAMETER THAT CHOOSES THE POLYGON.
%
% OUTPUT:
% [polygon_sides]: DEFINES THE COORDINATES OF THE VERTICES (REPEATING AT THE END
%                  THE FIRST POINT). THE VARIABLE "polygon_type" IS DEFINED ABOVE;
%                  "boundary_pts" IS A VECTOR HAVING AS COMPONENTS THE
%                  VERTICES OF THE COMPONENT, DESCRIBED COUNTERCLOCKWISE.
%                  POLYGON IS DEFINED INSIDE THE SQUARE [0,1] x [0,1].
%
%-------------------------------------------------------------------------------
% OBSERVE THAT THE FIRST VERTEX IS REPEATED (TO CLOSE THE POLYGON).
%-------------------------------------------------------------------------------

switch polygon_type
    
    case 0
        k_arc=64;
        k1=k_arc; % NUMBER OF SUBDIVISIONS 1ST ARC ECLIPSE.
        k2=k_arc; % NUMBER OF SUBDIVISIONS 2ND ARC ECLIPSE.
        
        a(1)=-pi/2; a(2)=pi;
        h1=(a(2)-a(1))/k1; v1=(a(1):h1:a(2))';
        polygon_sides_1=[0.5*(cos(v1)+1) 0.5*(sin(v1)+1)];
        
        a(3)=pi/2; a(4)=0;
        h2=(a(3)-a(4))/k1; v2=(a(3)-h2:-h2:a(4))';
        polygon_sides_2=[0.5*cos(v2) 0.5*sin(v2)];
        
        polygon_sides=[polygon_sides_1; polygon_sides_2];
        
        fprintf('\n \t [POLYGON]: ECLYPSE 2, %4.0f, %4.0f',...
            size(polygon_sides_1,1),size(polygon_sides,1));
        
        % plot(polygon_sides(:,1),polygon_sides(:,2),'r-');
        
    case 1
        fprintf('\n \t [POLYGON]: UNIT SQUARE [0,1]^2');
        k=1;
        polygon_sides=k*[0 0; 1 0; 1 1; 0 1; 0 0];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 2
        fprintf('\n \t [POLYGON]: CONVEX POLYGON');
        polygon_sides=[0.1 0; 0.7 0.2; 1 0.5; 0.75 0.85; 0.5 1; 0 0.25; 0.1 0];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 3
        fprintf('\n \t [POLYGON]: NON CONVEX POLYGON');
        polygon_sides=(1/4)*[1 2; 1 0; 3 2; 3 0; 4 2; 3 3; 3 0.85*4; 2 4; 0 3; 1 2];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 4
        fprintf('\n \t [POLYGON]: POLYGON');
        polygon_sides=(1/4)*[1 0; 3 2; 3 0; 4 2; 3 3; 3 4; 2 4; 0 3; 1 2; 1 0];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 5
        fprintf('\n \t [POLYGON]: POLYGON');
        polygon_sides=(1/4)*[0 0; 1 2; 2 0; 3 2; 4 0; 4 4; 3 2; 2 4; 1 2; 0 4; 0 0];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 6
        fprintf('\n \t [POLYGON]: UNIT SQUARE [-1,1]^2');
        polygon_sides=[-1 -1; 1 -1; 1 1; -1 1; -1 -1];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 7
        fprintf('\n \t [POLYGON]: UNIT TRIANGLE');
        polygon_sides=[0 0; 1 0; 1 1; 0 0];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 8
        k_arc=64;
        k1=k_arc; % NUMBER OF SUBDIVISIONS 1ST ARC ECLIPSE.
        k2=k_arc; % NUMBER OF SUBDIVISIONS 2ND ARC ECLIPSE.
        
        a(1)=-pi/2; a(2)=pi; a(3)=3*pi/2;
        h1=(a(2)-a(1))/k1; v1=(a(1):h1:a(2))';
        h2=(a(3)-a(2))/k2; v2=(a(2):h2:a(3))';
        polygon_sides_1=[cos(v1) sin(v1)];
        
        % BUILDING OPPOSITE CONCAVITY.
        C=[cos(v2)  sin(v2)];  % ROW VECTORS!
        L=size(C,1);
        A=[C(1,:)];            % ROW VECTORS!
        B=[C(L,:)];            % ROW VECTORS!
        v=B-A; v=v/norm(v,2);
        polygon_sides_2=[];
        for index=2:L
            C_LOC=C(index,:);
            u=C_LOC-A;
            u1=-u+2*(u*v')*v;
            polygon_sides_2=[polygon_sides_2; A+u1];
        end
        
        polygon_sides=[polygon_sides_1; polygon_sides_2];
        
        polygon_sides=0.5*ones(size(polygon_sides))+0.5*polygon_sides;
        fprintf('\n \t [POLYGON]: ECLYPSE, %4.0f, %4.0f',...
            size(polygon_sides_1,1),size(polygon_sides,1));
        
    case 9
        figure('position', get(0,'screensize'));
        axes('position', [0 0 1 1]);
        [x,y]=ginput(50);
        polygon_sides=[x y; x(1) y(1)];
        fprintf('\n \t [POLYGON]: MANUAL ENTRIES, %4.0f, %4.0f',...
            size(polygon_sides_1,1),size(polygon_sides,1));
        
    case 10
        
        k1=8; % NUMBER OF SUBDIVISIONS 1ST ARC CIRCLE.
        
        a(1)=0; a(2)=2*pi;
        h1=(a(2)-a(1))/k1; v1=(a(1):h1:a(2))';
        polygon_sides=[cos(v1) sin(v1)];
        
        polygon_sides=0.5*ones(size(polygon_sides))+0.5*polygon_sides;
        fprintf('\n \t [POLYGON]: DISK, %4.0f',size(polygon_sides,1));
    case 11
        polygon_sides=[1 0; 0.75 0.5; 0 1; 0 0.5; 1 0];
        fprintf('\n \t [POLYGON]: BOOMERANG, %4.0f',size(polygon_sides,1));
    case 12
        polygon_sides=[1 0; 0 1; -1 0; 0 -1; 1 0];
        fprintf('\n \t [POLYGON]: QUADRANGLE, %4.0f',size(polygon_sides,1));
    case 13
        fprintf('\n \t [POLYGON]: NON CONVEX POLYGON, HORSE LIKE');
        polygon_sides=(1/4)*[1 2; 1 0; 3 2; 3 0; 4 2; 3 3; 3 0.85*4; 4*0.6 4*0.95; 2 4; 0 3; 1 2];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
    case 14
        fprintf('\n \t [POLYGON]: NON CONVEX POLYGON, HAND');
        polygon_sides=[   0.49804687500000   0.23832417582418
            0.52148437500000   0.29601648351648
            0.55507812500000   0.35233516483516
            0.58398437500000   0.40315934065934
            0.58242187500000   0.42513736263736
            0.56757812500000   0.43475274725275
            0.54414062500000   0.41826923076923
            0.52226562500000   0.38942307692308
            0.49570312500000   0.38392857142857
            0.48710937500000   0.41552197802198
            0.48945312500000   0.47184065934066
            0.49882812500000   0.54601648351648
            0.50195312500000   0.62019230769231
            0.49257812500000   0.64766483516484
            0.47773437500000   0.63667582417582
            0.46601562500000   0.55975274725275
            0.45039062500000   0.47733516483516
            0.45273437500000   0.57760989010989
            0.45117187500000   0.65728021978022
            0.43085937500000   0.67788461538462
            0.41757812500000   0.65178571428571
            0.41679687500000   0.57074175824176
            0.41445312500000   0.48008241758242
            0.40195312500000   0.55700549450549
            0.38632812500000   0.63667582417582
            0.37304687500000   0.64766483516484
            0.36210937500000   0.62568681318681
            0.36992187500000   0.53640109890110
            0.37695312500000   0.44848901098901
            0.34492187500000   0.50068681318681
            0.32070312500000   0.54876373626374
            0.30195312500000   0.55151098901099
            0.29882812500000   0.53090659340659
            0.32460937500000   0.45398351648352
            0.34570312500000   0.40041208791209
            0.34960937500000   0.32074175824176
            0.35976562500000   0.23832417582418
            0.49804687500000   0.23832417582418];
        fprintf(' : %4.0f VERTICES',size(polygon_sides,1));
end








function z=fct2D(x,y,function_type,degree)

%-------------------------------------------------------------------------------
% FUNCTION DEFINITION.
%-------------------------------------------------------------------------------
%
% INPUT:
%
% [x,y]: EVALUATE THE FUNCTION IN THE POINTS (x,y). "x", "y" ARE COLUMN VECTORS.
%
% [function_type]: PARAMETER THAT CHOOSES THE FUNCTION.
%
% OUTPUT:
%
% [z]: VALUE OF THE FUNCTION IN THE POINTS (x,y). IT IS A COLUMN VECTOR.
%
%-------------------------------------------------------------------------------

if nargin < 4
    degree=20;
end

switch function_type
    case 0
        a=7.671018281063714e-01; b=3.547875152618030e-02;
        z= (1+a*(x-0.5) + b*(y-0.25)).^degree;
    case 1
        z= .75*exp(-((9*x-2).^2 + (9*y-2).^2)/4) + ...
            .75*exp(-((9*x+1).^2)/49 - (9*y+1)/10) + ...
            .5*exp(-((9*x-7).^2 + (9*y-3).^2)/4) - ...
            .2*exp(-(9*x-4).^2 - (9*y-7).^2);
    case 2
        z=( (x-0.5).^2 +(y-0.5).^2 ).^(1/2);
    case 3
        k=19;
        z=(x+y).^k;
    case 4
        loc_arg=(x-0.5).^2+(y-0.5).^2;
        z=exp(-loc_arg);
    case 5
        loc_arg=(x-0.5).^2+(y-0.5).^2;
        z=exp(-100*loc_arg);
    case 6
        z=cos(degree*(x+y));
    case 7
        z=ones(size(x));
    case 8
        z=exp(x+y);
    case 9
        den=1+16*(x.^2+y.^2);
        z=1./den;
    case 10
        exponent=3/2;
        z=(x.^2+y.^2).^(exponent);
end





%--------------------------------------------------------------------------
% Attached functions.
%--------------------------------------------------------------------------

function plot_spl(Sx,Sy,P,Z,Zin)

% ........................ plot preferences ...............................
plot_axis=0;

% ........................... plot: start .................................
clf;
hold on;

if plot_axis == 0
    axis off;
end

% ........................ plot: grids .........................

plot(Z(:,1),Z(:,2),'ro');
plot(Zin(:,1),Zin(:,2),'go');
plot(P(:,1),P(:,2),'k.');

% ....................... plot: parametric splines ........................

L=length(Sx);
for ii=1:L
    
    SxL=Sx(ii);
    SyL=Sy(ii);
    
    SxL_breaks=SxL.breaks;
    
    
    Nbreaks=length(SxL_breaks);
    
    for kk=1:Nbreaks-1
        
        t0=SxL_breaks(kk); t1=SxL_breaks(kk+1);
        
        N=10^3;
        tt=linspace(t0,t1,N);
        
        xx=ppval(SxL,tt);
        yy=ppval(SyL,tt);
        
        
        plot(xx,yy,'k-','LineWidth',1);
    end
    
end

hold off;


