
function demoS2_cub(domain_type)

%--------------------------------------------------------------------------
% Object:
% Demo on algebraic cubature rules some subdomains of the unit-sphere.
%--------------------------------------------------------------------------
% Dates:
% Written on 29/10/2020: M. Vianello;
%
% Modified on:
% 19/06/2023: A. Sommariva.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Domain to be considered. The variable "domain_type" can be:
%
% 1: 'sphere';
% 2: 'spherical-rectangle';
% 3: 'spherical-triangle';
% 4: 'spherical-polygon'.
%
% In case domain_type is 2, 3, 4, special examples are considered (for
% details, see the subroutine "define_domain" defined below).
%--------------------------------------------------------------------------

if nargin < 1, domain_type=2; end

%--------------------------------------------------------------------------
% Function to study. 
% The variable "gallery_str" can be 'easy' or 'hard'.
%--------------------------------------------------------------------------

gallery_str='easy'; example=8;

%--------------------------------------------------------------------------
% Degrees in numerical experiments: can be a vector.
%--------------------------------------------------------------------------

nV=2:2:10;


% ........................ Main code below ................................

comp_method=2; % 1. lsqnonneg, 2. LHDM (default: 2).

% ....... Apply settings to define domain, pointsets, and functions .......

% Domain
switch domain_type
    case 1, domain_example='sphere';
    case 2, domain_example='spherical-rectangle';
    case 3, domain_example='spherical-triangle';
    otherwise, domain_example='spherical-polygon';
end

% ........................ Main code below ................................

domain_struct=gallery_domains_S2(domain_example);
[g,gstr]=gallery_3D(gallery_str,example);

% ........ reference value ........

degR=max(nV)+20;

fprintf('\n \t * Computing Reference rule, ADE: %3.0f',degR);

tic; 
[XWR,dboxR]=define_cub_rule(domain_struct,degR); 
gXR=feval(g,XWR(:,1),XWR(:,2),XWR(:,3)); wR=XWR(:,4);
IR=wR'*gXR;
cpuR=toc;

fprintf('\n \t * Computing algebraic rules (full & compressed versions)');

for k=1:length(nV)
    
    n=nV(k);
    
  fprintf('\n \t -> ADE: %3.0f',n);


    % ........ full rule ........
    tic; [XW,dbox0]=define_cub_rule(domain_struct,n); cpu(k,1)=toc;
    gX=feval(g,XW(:,1),XW(:,2),XW(:,3)); w=XW(:,4);
    I(k)=w'*gX;
    AE(k)=abs(I(k)-IR); RE(k)=AE(k)./abs(IR);
    cardX(k)=size(XW,1);
    
    % ........ compressed rule ........
    X=XW(:,1:3); u=w;
    tic; [XC,wc,momerr(k),dbox1]=dCATCH(n,X,w,dboxR',[],(n+1)^2,comp_method); cpu(k,2)=toc;
    gXC=feval(g,XC(:,1),XC(:,2),XC(:,3));
    IC(k)=wc'*gXC;
    AEC(k)=abs(IC(k)-IR); REC(k)=AEC(k)./abs(IR);
    cardXC(k)=size(XC,1);
    
end



fprintf('\n \t * Statistics');

% ..... statistics ....
fprintf('\n \t ..............................................................................................');
fprintf('\n \t Domain: '); disp(domain_example);
fprintf('\n \t Function: '); disp(gstr);
fprintf('\n \t ..............................................................................................');
fprintf('\n \t |  n  | card X | cardXC |   cpuf   |   cpuc   |  momerr  |');
fprintf('\n \t ..............................................................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %6.0f | %6.0f | %1.2e | %1.2e | %1.2e |',...
        nV(k),cardX(k),cardXC(k),cpu(k,1),cpu(k,2),momerr(k));
end
fprintf('\n \t ..............................................................................................');
fprintf('\n \t | %3.0f | %6.0f | %6.0f | %1.2e | %1.2e |',...
    degR,length(XWR(:,1)),0,cpuR,0);
fprintf('\n \t ..............................................................................................');
fprintf('\n \n');

fprintf('\n \t ..............................................................................................');
fprintf('\n \t |  n  |          If           |          Ic           |   REf   |   REc   |  ');
fprintf('\n \t ..............................................................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %1.15e | %1.15e | %1.1e | %1.1e |',...
        nV(k),I(k),IC(k),RE(k),REC(k));
end
fprintf('\n \t ..............................................................................................');
fprintf('\n \t | %3.0f | %1.15e | %1.15e | %1.1e | %1.1e |',...
    degR,IR,IR,0,0);
fprintf('\n \t ..............................................................................................');
fprintf('\n \n');


% % ..... plots .....

fprintf('\n \t * Plots \n \n');

clear_figure(1)
figure(1);
title_str='Cubature pointset';
plot_S2(domain_struct,X,XC,title_str);
























