
function demo_spherical_polygons_QMC

% Modified: November 12, 2024.

example=1;
function_type=7;
deg_rule=3;

N_QMC=100000;
pointset_type='H';

[vertices,domain_str]=define_domain(example);
[f,fs]=define_function(function_type);

vertices=vertices./(vecnorm(vertices'))';


atol=10^(-6);
rtol=atol;

fprintf('\n * adaptive_cub_sphpgon');
[I,Ierr,flag,Ihigh,iters,tri_vertices,tri_conn_list,L1_vertices]=...
    adaptive_cub_sphpgon(vertices,f,atol,rtol);

fprintf('\n * cub_sphpgon');
[XW,tri_vertices,tri_conn_list]=cub_sphpgon(deg_rule,vertices);
fX=feval(f,XW(:,1),XW(:,2),XW(:,3));
Inum=(XW(:,4))'*fX;

fprintf('\n * cub_sphpgon_QMC');
XW_QMC=cub_sphpgon_QMC(N_QMC,vertices,pointset_type);
fXW_QMC=feval(f,XW_QMC(:,1),XW_QMC(:,2),XW_QMC(:,3));
IQMC=(XW_QMC(:,4))'*fXW_QMC;

fprintf('\n \t .....................');
fprintf('\n \t '); disp(fs);
fprintf('\n \t '); disp(domain_str);
fprintf('\n \n \t ..... geometry .....');
fprintf('\n \t vertices : %8.0f',size(tri_vertices,1)); 
fprintf('\n \t triangle : %8.0f',size(tri_conn_list,1)); 
fprintf('\n \n \t ..... adaptive .....');
fprintf('\n \t iters    : %8.0f',iters); 
fprintf('\n \t triangles: %8.0f',length(L1_vertices)); 
fprintf('\n \n \t ...... rule .......');
fprintf('\n \t pts      : %8.0f',size(XW,1)); 
fprintf('\n \n \t ..... errors ......');
fprintf('\n \t Ia: %-1.15e',I);
fprintf('\n \t In: %-1.15e',Inum);
fprintf('\n \t IQ: %-1.15e',IQMC);
fprintf('\n \n \t ... rule vs adp ...');
fprintf('\n \t AE: %-1.15e',I-Inum);
fprintf('\n \t In: %-1.15e',(I-Inum)/I);
fprintf('\n \n \t ... QMC vs adp ....');
fprintf('\n \t AE: %-1.15e',I-IQMC);
fprintf('\n \t In: %-1.15e',(I-IQMC)/I);
fprintf('\n \t ..................... \n \n');


clf
figure(1)
for k=1:size(tri_conn_list,1)
    con_L=tri_conn_list(k,:);
    tri_L=tri_vertices(con_L',:);
    tri_L(end+1,:)=tri_L(1,:);
    plot3(tri_L(:,1),tri_L(:,2),tri_L(:,3));
    hold on;
end
hold off;

figure(2)
for k=1:length(L1_vertices)
    tri_L=L1_vertices{k};
    tri_L(end+1,:)=tri_L(1,:);
    plot3(tri_L(:,1),tri_L(:,2),tri_L(:,3));
    hold on;
end
hold off;

figure(3)
plot3(vertices(:,1),vertices(:,2),vertices(:,3),'mo',...
    'MarkerFaceColor','b', ...
    'MarkerSize',2);
hold on;
plot3(XW(:,1),XW(:,2),XW(:,3),'mo',...
    'MarkerFaceColor','m', ...
    'MarkerSize',2);
hold off;






function [vertices,domain_str]=define_domain(example)

%--------------------------------------------------------------------------
% Object:
% This routine, defines the domain to analyse.
%--------------------------------------------------------------------------
% Input:
% example: determines the spherical triangle domain to be analysed. From 0
%   to 7, going from large domains to tiny ones (the higher "example", the
%   smaller the domain.
%--------------------------------------------------------------------------
% Output:
% vertices: it is a matrix whose rows are the vertices of the spherical
%           triangles.
% domain_str: string that stores the geometry of the region.
%--------------------------------------------------------------------------


switch example

    case 0
        % Define spherical polygon.
        factor=0.5;
        t=linspace(0,2*pi,10); t=(t(1:end-1))';
        XV=factor*cos(t); YV=factor*sin(t); ZV=sqrt(1-XV.^2-YV.^2);
        % Note: spherical polygon (first vertices not equal to last).
        vertices=[XV YV ZV];

        domain_str='Polygonal cap at North-Pole';

    case 1

        spec_settings=32;
        th=linspace(0,2*pi,spec_settings+1);
        %th=(th(1:end-1))';
        polygon_sides=[cos(th').*(1-cos(th')) sin(th').*(1-cos(th'))];
        polygon_sides=polygon_sides(1:end-1,:);
        XV=polygon_sides(:,1); XV=XV/2.1;
        YV=polygon_sides(:,2); YV=YV/2.1;
        ZV=sqrt(1-XV.^2-YV.^2);
        vertices=[XV YV ZV];

        domain_str='Polygonal cardioid at North-Pole';

    case 2

        % Australia island vertices in degrees.
        % Column 1: Longitude. Column 2: Latitude.
        vertices_degs=100*[1.462300280000000  -0.386970830000000
            1.451162780000000  -0.381483330000000
            1.396121110000000  -0.361567220000000
            1.385678610000000  -0.348266670000000
            1.374529440000000  -0.349083330000000
            1.359337220000000  -0.345341940000000
            1.340596110000000  -0.329141390000000
            1.235266670000000  -0.339379720000000
            1.150004170000000  -0.335283610000000
            1.131528610000000  -0.261491940000000
            1.142312220000000  -0.263158890000000
            1.168041670000000  -0.205237780000000
            1.221729170000000  -0.172633330000000
            1.239712780000000  -0.168258330000000
            1.243913060000000  -0.163392220000000
            1.252495830000000  -0.155816670000000
            1.259183610000000  -0.146770830000000
            1.268749720000000  -0.137471390000000
            1.283721390000000  -0.154966940000000
            1.292900000000000  -0.148620830000000
            1.303179720000000  -0.133733330000000
            1.308121390000000  -0.124083330000000
            1.317678890000000  -0.113258330000000
            1.347483890000000  -0.119512500000000
            1.360166940000000  -0.124987780000000
            1.359095560000000  -0.132850000000000
            1.395983610000000  -0.175387780000000
            1.415862220000000  -0.125583610000000
            1.437729170000000  -0.143991940000000
            1.462604440000000  -0.188633890000000
            1.500416670000000  -0.221254170000000
            1.531929440000000  -0.259316940000000
            1.525429440000000  -0.324442220000000
            1.499796110000000  -0.375050560000000
            1.462300280000000  -0.386970830000000];

        % Australia island vertices in radians.
        longitudes=deg2rad(vertices_degs(:,1));
        latitudes=deg2rad(vertices_degs(:,2));

        % Australia island vertices in cartesian coordinates.
        [XX,YY,ZZ] = sph2cart(longitudes,latitudes,1);
        vertices=[XX YY ZZ];

        domain_str='Australia as island (33 vertices, via data)';

    case 3


        % Australia island (no Tasmania), vertices in degrees.
        % Column 1: Longitude. Column 2: Latitude.
        % Planar domain in polyshape form.
        australia_pshape=coastline_australia(0);
        vertices_degs=australia_pshape.Vertices;

        % Australia & Tasmania vertices in radians.
        longitudes=deg2rad(vertices_degs(:,1));
        latitudes=deg2rad(vertices_degs(:,2));

        % Australia & Tasmania vertices in cartesian coordinates.
        [XX,YY,ZZ] = sph2cart(longitudes,latitudes,1);
        vertices=[XX YY ZZ];

        domain_str='Australia as island (170 vertices, via polyshape)';

    case 4
        vertices = [0.4789    0.3966   -0.7832
            0.2813   -0.5443   -0.7904
            0.1032   -0.9815    0.1615
            0.3007    0.7574    0.5796];

        domain_str='Quadrangle';

    case 5
        vertices =[-0.4886   -0.5505    0.6769
            -0.3390   -0.4644    0.8182
            -0.3114   -0.6887   -0.6547
            -0.5290   -0.3708   -0.7633
            0.5768    0.7093   -0.4053
            -0.9359   -0.2140   -0.2799
            -0.6399   -0.7358    0.2216];

        domain_str='Eptagon';

end








function [f,fs]=define_function(function_type,parms)

%--------------------------------------------------------------------------
% Object:
% This routine, defines a function "f" to approximate and a string "fs" for
% possible messages to the user.
%--------------------------------------------------------------------------
% Input:
% function_type: determines the function to study.
% The first five functions has been used in the paper mentioned below.
%
% case 1, f=@(x,y,z) 1+x+y.^2+x.^2.*y+x.^4+y.^5+x.^2.*y.^2.*z.^2;
% case 2, f=@(x,y,z) cos(10*(x+y+z));
% case 3  x0=0; y0=0; z0=1;
%         g=@(x,y,z) (x-x0).^2 + (y-y0).^2 + (z-z0).^2;
%         f=@(x,y,z) exp( - g(x,y,z) );
% case 4  centroid=sum(vertices,1); centroid=centroid/norm(centroid);
%         x0=centroid(1); y0=centroid(2); z0=centroid(3);
%         g=@(x,y,z) (x-x0).^2 + (y-y0).^2 + (z-z0).^2;
%         f=@(x,y,z) exp( - g(x,y,z) );
% case 5  x0=0; y0=0; z0=1;
%         f=@(x,y,z) ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2);
% case 6  centroid=sum(vertices,1); centroid=centroid/norm(centroid);
%         x0=centroid(1); y0=centroid(2); z0=centroid(3);
%         f=@(x,y,z) ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2);
%--------------------------------------------------------------------------
% Output:
% f: defines a function "f" to approximate;
% fs: string with the content of the function "f".
%--------------------------------------------------------------------------
% Reference paper:
% A. Sommariva and M. Vianello
% Numerical hyperinterpolation overspherical triangles
%--------------------------------------------------------------------------

switch function_type
    case 1 % Fornberg
        f=@(x,y,z) 1+x+y.^2+x.^2.*y+x.^4+y.^5+x.^2.*y.^2.*z.^2;
        fs='1+x+y.^2+x.^2.*y+x.^4+y.^5+x.^2.*y.^2.*z.^2';

    case 2
        f=@(x,y,z) cos(10*(x+y+z));
        fs='cos(10*(x+y+z))';

    case 3 % exp and - square north pole distance
        x0=0; y0=0; z0=1;
        g=@(x,y,z) (x-x0).^2 + (y-y0).^2 + (z-z0).^2;
        f=@(x,y,z) exp( - g(x,y,z) );
        fs='exp(-g(x,y,z)), g=@(x,y,z) (x-x0).^2+(y-y0).^2+(z-z0).^2, x0=0;y0=0;z0=1;';

    case 4 % exp and - square centroid distance
        if nargin < 1
            vertices=[0 0 1; 1 0 0; 0 1 0];
        else
            vertices=parms;
        end

        % ... centroid computation ...
        ok_vertices=find( isnan(vertices(:,1))== 0 & isnan(vertices(:,2)) == 0 );
        vertices=vertices(ok_vertices,:);
        centroid=sum(vertices,1); centroid=centroid/norm(centroid);
        x0=centroid(1); y0=centroid(2); z0=centroid(3);

        % ... function ...
        g=@(x,y,z) (x-x0).^2 + (y-y0).^2 + (z-z0).^2;
        f=@(x,y,z) exp( - g(x,y,z) );
        fs='exp(-g(x,y,z)), g=@(x,y,z) (x-x0).^2+(y-y0).^2+(z-z0).^2';

        % ... output string ...
        x0str=num2str(x0,'%1.3e');
        y0str=num2str(y0,'%1.3e');
        z0str=num2str(z0,'%1.3e');
        fs='exp(-g(x,y,z)), g=@(x,y,z) (x-x0).^2+(y-y0).^2+(z-z0).^2';
        fs=strcat(fs,'  centroid=(',x0str,',',y0str,',',z0str,')');

    case 5 % north pole distance like
        x0=0; y0=0; z0=1;
        f=@(x,y,z) ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2);
        fs='((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2), x0=0; y0=0; z0=1;';

    case 6 % centroid distance like
        if nargin < 1
            vertices=[0 0 1; 1 0 0; 0 1 0];
        else
            vertices=parms;
        end

        % ... centroid computation ...
        ok_vertices=find( isnan(vertices(:,1))== 0 & ...
            isnan(vertices(:,2)) == 0 );
        vertices=vertices(ok_vertices,:);
        centroid=sum(vertices,1); centroid=centroid/norm(centroid);
        x0=centroid(1); y0=centroid(2); z0=centroid(3);

        % ... function ...
        f=@(x,y,z) ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2);

        % ... output string ...
        x0str=num2str(x0,'%1.3e');
        y0str=num2str(y0,'%1.3e');
        z0str=num2str(z0,'%1.3e');
        fs='((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2), ';
        fs=strcat(fs,'  centroid=(',x0str,',',y0str,',',z0str,')');

    case 7

        x0=-0.4886;   y0=-0.5505;    z0=0.6769;
        f=@(x,y,z) exp(-((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2))+...
            sin(100*x+5*y-0.7*z.^2);
        fs='exp(-((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(5/2))+...';
end  % end: "define_function"


