
function [Ntri,Nrule,NruleC,momerr,cpus]=demo_cub_continents(...
    continent_parm,deg,doplot)

%--------------------------------------------------------------------------
% Object:
%
% Demo of cubature over spherical polygons, in which the domains are
% coarse approximation of continents.
%--------------------------------------------------------------------------
% Input:
%
% continent_parm   : 1. South/North America, 2. Africa, 3 or 4. Australia.
% deg              : Algebraic degree of exactness.
% doplot           : 0: no plots, 1: plots.
%--------------------------------------------------------------------------
% Output:
%
% Ntri             : Number of spherical triangles.
% Nrule            : Cardinality of the rule.
% NruleC           : Cardinality of the compressed rule.
% momerr           : Moment matching.
% cpus             : Array of cputimes
%                    cpus(1): polyshape + triangulation
%                    cpus(2): basic rule computation
%                    cpus(3): compressed rule computation (from basic rule)
%--------------------------------------------------------------------------
% Information:
% Authors: A. Sommariva and M. Vianello.
% Update: January 12, 2022.
%--------------------------------------------------------------------------
%% Copyright (C) 2021
%% Alvise Sommariva, Marco Vianello.
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Authors:
%% Alvise Sommariva, Marco Vianello.
%%
%% Date: JULY 13, 2021-
%% Significant update: January 12, 2022.
%--------------------------------------------------------------------------

% clf; clear all; close all;

if nargin < 1
    % 1. South/North America, 2. Africa, 3 or 4. Australia, 5: Africa cont.
    continent_parm=5;
end

if nargin < 2
    % Algebraic degree of exactness of the cubature rule.
    deg=10;
end

if nargin < 3
    doplot=1;        % 0: no plots, 1: plots.
end
write_stats=0;   % 0: no stats, 1: stats.







%---------------------------- Main code below------------------------------

% In vertices we store, respectively, Longitude and Latitute, in degrees.
switch continent_parm
    case 1
        str='South-America and North-America';
        continent_pshape=coastline_america;

    case 2
        str='Africa';
        continent_pshape=coastline_africa;

    case 3
        str='Australia island';
        continent_pshape=coastline_australia(0);

    case 4

        % As an example, we add directly the data of Longitudes and
        % Latitutes in degrees, so that one can easily adapt the code to
        % other domains.

        str='Australia coarse (33 vertices)';
        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];

        continent_pshape=polyshape(vertices_degs);

    case 5


        vertices_degs=[
    32.478361000000000 29.940417000000000
    32.625388999999998 28.976638999999999
    33.586221999999999 27.847500000000000
    35.138750000000002 24.513332999999999
    35.793749999999996 23.909222000000000
    35.476222000000000 23.934221999999998
    35.667888999999995 22.946610999999997
    36.895416999999995 22.063333000000000
    37.318750000000001 21.064999999999998
    37.111193999999998 21.214167000000000
    37.412056000000000 18.875000000000000
    38.557055999999996 18.060110999999999
    39.712027999999997 15.094999999999999
    39.875000000000000 15.507888999999999
    40.209972000000000 14.948722000000000
    41.189971999999997 14.623721999999999
    43.327027999999999 12.475778000000000
    43.376193999999998 11.988360999999999
    42.527471999999996 11.522888999999999
    43.156610999999998 11.612943999999999
    44.575806000000000 10.375389000000000
    45.805833000000000 10.867889000000000
    46.440805999999995 10.683750000000000
    47.412499999999994 11.179582999999999
    49.426666999999995 11.337888999999999
    50.797610999999996 11.986110999999999
    51.291249999999998 11.833306000000000
    51.022055999999999 10.420805999999999
    51.415416999999998 10.443332999999999
    50.902055999999995 10.314943999999999
    50.830416999999997 9.419167000000000
    47.965443999999998 4.477500000000000
    46.011666999999996 2.429528000000000
    43.559166999999995 0.714583000000000
    41.314971999999997 -1.961250000000000
    40.774555999999997 -1.934167000000000
    40.987916999999996 -2.264194000000000
    40.186278000000001 -2.735833000000000
    39.201943999999997 -4.680861000000000
    38.772888999999999 -6.052528000000000
    39.551249999999996 -7.017500000000000
    39.217110999999996 -7.850833000000000
    39.447055999999996 -7.813361000000000
    39.255389000000001 -8.292527999999999
    39.704443999999995 -10.031832999999999
    40.642055999999997 -10.684277999999999
    40.332082999999997 -11.313333000000000
    40.642083000000000 -12.746694000000000
    40.394556000000001 -12.928388999999999
    40.842110999999996 -14.832528000000000
    39.856639000000001 -16.438806000000000
    36.854971999999997 -17.878750000000000
    36.312500000000000 -18.877110999999999
    34.884999999999998 -19.852944000000001
    34.631222000000001 -19.650832999999999
    34.661221999999995 -20.555056000000000
    35.546222000000000 -22.171721999999999
    35.502055999999996 -24.107527999999999
    32.498666999999998 -25.960055999999998
    32.834055999999997 -26.292916999999999
    32.964388999999997 -26.084111000000000
    32.385388999999996 -28.549194000000000
    31.274555999999997 -29.454193999999998
    30.009166999999998 -31.297110999999997
    27.011638999999999 -33.565416999999997
    24.839971999999999 -34.206277999999998
    22.555833000000000 -33.978749999999998
    20.019971999999999 -34.830444000000000
    18.790806000000000 -34.080444000000000
    18.399943999999998 -34.302110999999996
    17.847888999999999 -32.820028000000001
    18.327082999999998 -32.530028000000001
    18.200443999999997 -31.682499999999997
    14.946249999999999 -26.317527999999999
    14.502027999999999 -22.542610999999997
    11.799500000000000 -17.979194000000000
    11.744555999999999 -15.835832999999999
    12.532888999999999 -13.418360999999999
    13.657055999999999 -12.235028000000000
    13.862833000000000 -11.008360999999999
    12.994582999999999 -9.087527999999999
    13.392111000000000 -8.388332999999999
    12.272917000000000 -6.140000000000000
    12.672943999999999 -6.010860999999999
    12.167916999999999 -5.727556000000000
    11.727860999999999 -4.450861000000000
    9.619583000000000 -2.408361000000000
    8.700388999999999 -0.668361000000000
    9.130417000000000 -0.714194000000000
    9.351638999999999 0.357917000000000
    9.501666999999999 0.050389000000000
    9.987083000000000 0.172500000000000
    9.309583000000000 0.512500000000000
    9.663750000000000 0.452500000000000
    9.695971999999999 1.078750000000000
    9.343722000000000 1.175833000000000
    9.975417000000000 3.079111000000000
    9.542916999999999 3.810833000000000
    9.833693999999999 3.809167000000000
    9.727943999999999 4.097500000000000
    9.214167000000000 3.941194000000000
    8.830389000000000 4.746667000000000
    8.554971999999999 4.482056000000000
    8.235778000000000 4.928750000000000
    8.278278000000000 4.536222000000000
    7.340000000000000 4.443722000000000
    7.075806000000000 4.737916999999999
    7.005777999999999 4.372917000000000
    6.779139000000000 4.667860999999999
    6.867471999999999 4.359583000000000
    6.723332999999999 4.612916999999999
    6.705028000000000 4.338722000000000
    5.983333000000000 4.310417000000000
    5.362083000000000 5.163306000000000
    5.680444000000000 5.540806000000000
    5.361667000000000 5.390416999999999
    5.176250000000000 5.568306000000000
    5.542860999999999 5.619193999999999
    5.080388999999999 5.718332999999999
    5.367083000000000 5.970000000000000
    5.064166999999999 5.764527999999999
    4.430833000000000 6.350389000000000
    3.401250000000000 6.395805999999999
    3.821639000000000 6.621194000000000
    1.343361000000000 6.161167000000000
    -1.981694000000005 4.754611000000000
    -4.016999999999996 5.258389000000000
    -6.124167000000000 4.997027999999999
    -7.480806000000030 4.360361000000000
    -8.717500000000030 4.807083000000000
    -11.230750000000000 6.789499999999999
    -12.521277999999995 7.384111000000000
    -12.456250000000011 7.766639000000000
    -13.293778000000032 8.421638999999999
    -12.837083000000007 8.560832999999999
    -13.174611000000027 8.531666999999999
    -13.303944000000001 9.334861000000000
    -13.723749999999995 9.506639000000000
    -13.537056000000007 9.796611000000000
    -14.455833000000041 10.208722000000000
    -14.697083000000021 11.076639000000000
    -15.078389000000016 10.918721999999999
    -15.004638999999997 11.183306000000000
    -15.235000000000014 10.997916999999999
    -15.504611000000011 11.328332999999999
    -15.033749999999998 11.643277999999999
    -15.509583000000021 11.789166999999999
    -14.939166999999998 11.988721999999999
    -15.955833000000041 11.734556000000000
    -15.853778000000034 11.992471999999999
    -16.352082999999993 12.096667000000000
    -16.014194000000032 12.345388999999999
    -16.776250000000005 12.406639000000000
    -16.220472000000029 12.609971999999999
    -16.763778000000002 12.571610999999999
    -16.717943999999989 13.459999999999999
    -16.355028000000004 13.249582999999999
    -15.574611000000004 13.515806000000000
    -16.506666999999993 13.361250000000000
    -16.513778000000002 13.879166999999999
    -17.533778000000041 14.741667000000000
    -16.577056000000027 15.697471999999999
    -16.043749999999989 17.747499999999999
    -16.547917000000041 19.379999999999999
    -16.197917000000018 20.224971999999998
    -16.912528000000009 21.161193999999998
    -17.107222000000036 20.837527999999999
    -16.983000000000004 21.798138999999999
    -15.704583000000014 23.972500000000000
    -15.915416999999991 23.803305999999999
    -14.903777999999988 24.685860999999999
    -14.465417000000002 26.196666999999998
    -13.651250000000005 26.659193999999999
    -12.915861000000007 27.956249999999997
    -11.525000000000034 28.298749999999998
    -10.217943999999989 29.309138999999998
    -9.662944000000039 30.071638999999998
    -9.846222000000012 31.397499999999997
    -9.287111000000039 32.545833000000002
    -6.812500000000000 34.025388999999997
    -5.934639000000004 35.790777999999996
    -5.470861000000014 35.917110999999998
    -4.688389000000029 35.206277999999998
    -2.979139000000032 35.441972000000000
    -2.830083000000002 35.108666999999997
    -1.957500000000039 35.075389000000001
    -0.937528000000043 35.727916999999998
    0.000000000000000 35.842889000000000
    1.084139000000000 36.493749999999999
    3.895861000000000 36.925388999999996
    5.344139000000000 36.641193999999999
    6.415806000000000 37.086250000000000
    7.899972000000000 36.842917000000000
    9.744999999999999 37.349556000000000
    10.431611000000000 36.712888999999997
    11.039138999999999 37.089582999999998
    10.467832999999999 36.135777999999995
    11.162917000000000 35.215832999999996
    10.017861000000000 34.187528000000000
    10.214138999999999 33.797083000000001
    12.299139000000000 32.837888999999997
    13.388306000000000 32.897860999999999
    15.205027999999999 32.382055999999999
    15.741667000000000 31.397888999999999
    17.853306000000000 30.925388999999999
    18.989971999999998 30.274583000000000
    20.057082999999999 30.850860999999998
    20.099999999999998 32.182055999999996
    21.708306000000000 32.941221999999996
    23.107471999999998 32.635388999999996
    23.319944000000000 32.152889000000002
    24.979832999999999 31.971055999999997
    25.200832999999999 31.523638999999999
    26.982471999999998 31.442055999999997
    29.064943999999997 30.822056000000000
    31.028305999999997 31.600389000000000
    31.919166999999998 31.531222000000000
    32.478361000000000 29.940417000000000];

        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];


        continent_pshape=polyshape(vertices_degs);
        domain_str='Africa (228 vertices, via data)';



end



% fprintf(2,'\n \t Continent: '); disp(str);

% fprintf(2,'\n \t Determine triangulation: ');
tic;
% Make polyshape and triangulation
tri = triangulation(continent_pshape);
vertices=tri.Points;
cpus(1)=toc;




% Determine full cubature rule
% fprintf(2,'\n \t Determine full cubature rule: ');
tic;
xyzw=[];
for k=1:size(tri,1)

    % Local triangle
    triL=tri(k,:);

    % Vertices from degrees to radians, to cartesian coordinates.
    P1deg=vertices(triL(1),:); P1rad=deg2rad(P1deg);
    [P1x,P1y,P1z] = sph2cart(P1rad(1),P1rad(2),1); P1=[P1x P1y P1z]';

    P2deg=vertices(triL(2),:); P2rad=deg2rad(P2deg);
    [P2x,P2y,P2z] = sph2cart(P2rad(1),P2rad(2),1); P2=[P2x P2y P2z]';

    P3deg=vertices(triL(3),:); P3rad=deg2rad(P3deg);
    [P3x,P3y,P3z] = sph2cart(P3rad(1),P3rad(2),1); P3=[P3x P3y P3z]';

    % Local rule.
    [xyzwL] = cub_sphtri(deg,P1,P2,P3,0);

    % Add to storage.
    xyzw=[xyzw; xyzwL];

end
cpus(2)=toc;




% Determine Caratheodory-Tchakaloff compression.
% fprintf(2,'\n \t Determine Caratheodory-Tchakaloff compression: ');
tic;
[nodes,w,momerr,dbox] = dCATCHsph(deg,xyzw(:,1:3),xyzw(:,4));
cpus(3)=toc;

Ntri=size(tri);  Nrule=size(xyzw,1); NruleC=size(nodes,1);

if write_stats
    fprintf('\n \t .................... Stats ...............................')
    fprintf('\n \t ADE: %3.0f',deg);
    fprintf('\n \t CARD TRI: %3.0f',Ntri);
    fprintf('\n \t CARD F  : %3.0f',Nrule);
    fprintf('\n \t CARD C  : %3.0f',NruleC);
    fprintf('\n \t MOM.ERR.: %1.2e',momerr);
    fprintf('\n \t CPUTIMES: %1.2e',sum(cpus));
    fprintf('\n \n')
end




% ............................ 2D plot ....................................

% Figure 1. Plot triangulation
if doplot


    % ................. Figure 1 .................
    % 1. Plot continent.
    figure(1)
    hold on;
    plot(continent_pshape,'EdgeColor','black','LineWidth',4);
    triplot(tri);


    % 2. Plot nodes (full rule)
    XX=xyzw(:,1); YY=xyzw(:,2); ZZ=xyzw(:,3);
    [TH,PHI,R] = cart2sph(XX,YY,ZZ);
    THdeg=rad2deg(TH);
    PHIdeg=rad2deg(PHI);
    plot(THdeg,PHIdeg,'g.');


    % % 3. Plot nodes (compressed rule)
    XXc=nodes(:,1); YYc=nodes(:,2); ZZc=nodes(:,3);
    [THc,PHIc,R] = cart2sph(XXc,YYc,ZZc);
    THdegc=rad2deg(THc);
    PHIdegc=rad2deg(PHIc);
    plot(THdegc,PHIdegc,'mo','MarkerEdgeColor','k',...
        'MarkerFaceColor','k',...
        'MarkerSize',8)

    % axis off

    % 4. Add title.
    title_str=strcat('Cubature nodes projection, ADE= ',num2str(deg));
    title(title_str);

    hold off;






    % ............................ 3D plot ................................



    % ................. Figure 2 .................
    % * Plot Full set

    figure(2);

    % 1. Plot set.
    hold on; plot_s2('sphere',[],xyzw); hold on;

    % 2. Add title.
    title_str=strcat('Cubature nodes: full set, ADE= ',num2str(deg));
    title(title_str);

    hold off;




    % ................. Figure 3 .................
    % * Plot Compressed set and triangulation

    figure(3)
    hold on;

    % 1. Plot sphere.
    [X,Y,Z] = sphere(20);
    R=max(sqrt(XX.^2+YY.^2+ZZ.^2));
    plot3(R*X,R*Y,R*Z,'MarkerFaceColor','red');

    % 2. Plot spherical polygon triangulation.
    for k=1:size(tri,1)

        triL=tri(k,:); % local triangle
        % from degrees to radians, to cartesian

        P1deg=vertices(triL(1),:); P1rad=deg2rad(P1deg);
        [P1x,P1y,P1z] = sph2cart(P1rad(1),P1rad(2),1); P1=[P1x P1y P1z];

        P2deg=vertices(triL(2),:); P2rad=deg2rad(P2deg);
        [P2x,P2y,P2z] = sph2cart(P2rad(1),P2rad(2),1); P2=[P2x P2y P2z];

        P3deg=vertices(triL(3),:); P3rad=deg2rad(P3deg);
        [P3x,P3y,P3z] = sph2cart(P3rad(1),P3rad(2),1); P3=[P3x P3y P3z];

        plot_side(P1,P2); plot_side(P2,P3); plot_side(P3,P1);

        hold on;
    end


    % 3. Plot compressed set.
    plot_s2('sphere',[],[],nodes);


    % 4. Add title.
    title_str=...
        strcat('Triangulation and Caratheodory-Tchakaloff nodes, ADE= ',...
        num2str(deg));
    title(title_str);
    hold off;
end










function plot_side(A,B)

%--------------------------------------------------------------------------
% Object:
% Plot "geodesic line" between A and B.
%--------------------------------------------------------------------------
% Input:
% A,B: row vectors of two points of the sphere to be joined.
%--------------------------------------------------------------------------

R=norm(A);
t=linspace(0,1,10);
C=[];
for k=1:length(t)
    CC=A+t(k)*(B-A); CC=CC/norm(CC); C=[C; R*CC];
end

mycolor='black';
plot3(C(:,1),C(:,2),C(:,3),mycolor,'LineWidth',1);

