
function demo_cub_sphpoly

clf; clear all;

addpath(genpath('../EXTERNAL_ROUTINES/'));

example=1;
N=1000;
pointset_type='H';

f=@(x,y,z) x.^2+7*y.^3+0.5*z.^4;
f=@(x,y,z) 1+1*x+0*y+0*z;

vertices = sph_domain(example);

xyzwR=cub_sphpgon(15,vertices);
xyzw=cub_sphpoly_QMC(N,vertices,pointset_type);

xR=xyzwR(:,4); yR=xyzwR(:,2); zR=xyzwR(:,3); wR=xyzwR(:,4);
fR=feval(f,xR,yR,zR);
IR=wR'*fR;

xQ=xyzw(:,1); yQ=xyzw(:,2); zQ=xyzw(:,3); wQ=xyzw(:,4);
fQ=feval(f,xQ,yQ,zQ);
IQ=wQ'*fQ;

fprintf('\n \t %1.15e',IR);
fprintf('\n \t %1.15e',IQ);

fprintf('\n')



% plot_S2('sphere',xyzw(:,1:3),xyzwR(:,1:3));

clf(1)
figure(1)
plot_S2('sphere',xyzwR(:,1:3));
hold on;
plot_sphpoly_bnd(vertices);
hold off;

clf(2)
figure(2)
plot_S2('sphere',xyzw(:,1:3));
hold on;
plot_sphpoly_bnd(vertices);
hold off;














function V = sph_domain(example)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Depending on the variable "example", the routine:
%
% a) determines the vertices of the N sides spherical polygon whose
% vertices are stored in the N x 3 matrix V.
%
% b) computes the area of the spherical domain.
%
% c) determines a set of uniform points on the spherical domain, starting
% from "N" points in the "minimal cap" containing the domain.
%--------------------------------------------------------------------------


% ........................ domain definition ..............................

switch example
    case 0
        V = [1 0 0; 0 1 0; 0 0 1];


    case 1
        V = [0.4789    0.3966   -0.7832
            0.2813   -0.5443   -0.7904
            0.1032   -0.9815    0.1615
            0.3007    0.7574    0.5796];
    case 2
        V =[-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];
    case 3
        australia_pshape = coastline_australia(0);

        Vdeg = australia_pshape.Vertices;
        [Vx,Vy,Vz] = sph2cart(deg2rad(Vdeg(:,1)),deg2rad(Vdeg(:,2)),1);
        V = [Vx,Vy,Vz];
    case 4
        america_pshape = coastline_america;

        Vdeg = america_pshape.Vertices;
        [Vx,Vy,Vz] = sph2cart(deg2rad(Vdeg(:,1)),deg2rad(Vdeg(:,2)),1);
        V = [Vx,Vy,Vz];
    case 5
        africa_pshape = coastline_africa;

        polyAfrica = regions(africa_pshape);
        poly1 = polyAfrica(2);
        Vdeg = poly1.Vertices;
        [Vx,Vy,Vz] = sph2cart(deg2rad(Vdeg(:,1)),deg2rad(Vdeg(:,2)),1);
        V = [Vx,Vy,Vz];
end