
function demo_sphgon_area

% IMPORTANT: Requires Mapping toolbox.

example=1;

% Define domain
[vertices,domain_str]=define_domain(example);

% Convert from Cartesian to Longitude-Latitude
[lon,lat,h]=cart2sph(vertices(:,1),vertices(:,2),vertices(:,3));
lon=lon+2*pi*(lon < 0);

a = areaint(lat,lon,'radians')

% Sph routine
[XW,tri_vertices,tri_conn_list]=cub_sphpgon(1,vertices);

f=@(x,y,z) 1+0*x;
fX=feval(f,XW(:,1),XW(:,2),XW(:,3));
Inum=(XW(:,4))'*fX

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


