
function [IH,IL,flag,iters]=cub_polygon_adaptive(vertices,f,tol)

%--------------------------------------------------------------------------
% Input:
% vertices : it is a M x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex of the polygonal
%            region (counterclockwise order);
% f        : function to integrate over the spherical triangle defined by
%            vertices;
% tol      : absolute cubature error tolerance.
%--------------------------------------------------------------------------
% Output:
% I      : approximation of integral of "f" over the spherical triangle
%          defined by vertices;
% IH     : high order approximation of the integral;
% IL     : low order approximation of the integral;
%--------------------------------------------------------------------------
% Subroutines:
% 1. add_triangle_data
% 2. generate_triangle_sons
% 3. cub_tri
% 4. lyness_jespersen
%--------------------------------------------------------------------------
% Dates:
% Initial version: January 05, 2021, by A. Sommariva
%--------------------------------------------------------------------------

% ................. Troubleshooting and settings ..........................

if nargin < 3, tol=10^(-6); end

% max number of triangles in which the domain is partitioned
max_triangles=5000;

% ........ Procedure start up, data stored in structures ............

%..........................................................................
% Main strategy:
% First we triangulate the spherical polygon, providing the vertices of the
% triangulation and its connectivity list (description of each spherical
% triangle via its vertices).
% We compute low and high order integrals on each triangle and then keep
% erasing that one with worst error from the list, subdivide it in
% simplices, compute integrals in each of them and updating the list.
%..........................................................................

flag=0;

L1_vertices={}; L1_integrals_L=[]; L1_integrals_H=[]; L1_errors=[];

% If input is not a polyshape determines its polyshape "PG".
S=class(vertices);
if strcmp(S,'polyshape') == 0
    XV=vertices(:,1); YV=vertices(:,2); PG=polyshape(XV,YV);
else
    PG=vertices;
end

tri=triangulation(PG);
tri_vertices=tri.Points; tri_conn_list=tri.ConnectivityList;

for k=1:size(tri_conn_list,1)
    verticesL=tri_vertices((tri_conn_list(k,:))',:);
    [L1_vertices,L1_integrals_L,L1_integrals_H,L1_errors]=...
        add_triangle_data(verticesL,f,L1_vertices,L1_integrals_L,...
        L1_integrals_H,L1_errors);
end

IL=sum(L1_integrals_L); IH=sum(L1_integrals_H); AE=abs(IL-IH);

% ........................ successive refinement ..........................

iters=1;

while (AE >= tol) & (AE >= tol*IH)
    
    N=length(L1_integrals_L);
    
    % Too many triangles: exit with errors
    if N > max_triangles, flag=1; return; end
    
    % Determone triangle with worst error.
    [Ierr_max,kmax]=max(L1_errors);
    vertices_kmax=L1_vertices{kmax};
    
    % Erase "kmax" triangle from LIST 1
    k_ok=setdiff(1:N,kmax);
    
    if length(k_ok) > 0
        L1_vertices={L1_vertices{k_ok}};
        L1_integrals_L=L1_integrals_L(k_ok);
        L1_integrals_H=L1_integrals_H(k_ok);
        L1_errors=L1_errors(k_ok);
    else
        L1_vertices={}; L1_integrals_L=[]; L1_integrals_H=[]; L1_errors=[];
    end
    
    new_triangles_cell=generate_triangle_sons(vertices_kmax);
    
    % Replace worst triangle "data", with "those" of a triangulation.
    for j=1:4
        verticesL=new_triangles_cell{j};
        [L1_vertices,L1_integrals_L,L1_integrals_H,L1_errors]=...
            add_triangle_data(verticesL,f,L1_vertices,L1_integrals_L,...
            L1_integrals_H,L1_errors); 
    end
    
    IL=sum(L1_integrals_L); IH=sum(L1_integrals_H); AE=abs(IL-IH);
 
    iters=iters+1;
    
end









function [L1_vertices,L1_integrals_L,L1_integrals_H,L1_errors]=...
    add_triangle_data(vertices,f,L1_vertices,L1_integrals_L,...
    L1_integrals_H,L1_errors)


L1_vertices{end+1}=vertices;
[IL,IH]=cub_tri(f,vertices);
L1_integrals_L(end+1)=IL;
L1_integrals_H(end+1)=IH;
L1_errors(end+1)=abs(IL-IH);









function triangles_cell=generate_triangle_sons(vertices)

%--------------------------------------------------------------------------
% Input:
% vertices : it is a 3 x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex (counterclockwise);
% f        : function to integrate over the spherical triangle defined by
%            vertices;
%--------------------------------------------------------------------------
% Output:
% L2_data :
%--------------------------------------------------------------------------

OA=(vertices(1,:)); OB=(vertices(2,:)); OC=(vertices(3,:));

% ........................ compute midpoints ..............................

OAB_mid=(OA+OB)/2; OAC_mid=(OA+OC)/2;  OBC_mid=(OB+OC)/2;

% ........................ triangle data ..................................


triangles_cell{1}=[OA;OAB_mid;OAC_mid];
triangles_cell{2}=[OAB_mid; OB; OBC_mid];
triangles_cell{3}=[OBC_mid; OC; OAC_mid];
triangles_cell{4}=[OAB_mid; OBC_mid; OAC_mid];









function [IL,IH]=cub_tri(f,vertices)

%--------------------------------------------------------------------------
% Object:
% Computation of low order and higher order rules of integral of "f" on the
% triangle T determined by vertices.
%--------------------------------------------------------------------------
% Input:
% f: function to integrate;
% vertices: 3 x 3 matrix whose k-th row represents the cartesian
% coordinates of the k-th vertex.
%--------------------------------------------------------------------------
% Output:
% IL: low order approximation of the integral of "f" on the triangle T;
% IH: high order approximation of the integral of "f" on the triangle T.
%--------------------------------------------------------------------------

tri_area=polyarea(vertices(:,1),vertices(:,2));

% ........................... low order rule ..............................

xyw_bar=lyness_jespersen(9);

% nodes and weights
lambda=xyw_bar(:,1:3);
nodes=lambda*vertices;

% Low integral
w=tri_area*xyw_bar(:,4); 

fXY=feval(f,nodes(:,1),nodes(:,2)); 
IL=w'*fXY;

% ........................... high order rule ..............................

xyw_bar=lyness_jespersen(11);

% nodes and weights
lambda=xyw_bar(:,1:3);
nodes=lambda*vertices;

% Low integral
w=tri_area*xyw_bar(:,4); fXY=feval(f,nodes(:,1),nodes(:,2)); IH=w'*fXY;









function xyw_bar=lyness_jespersen(deg)

switch deg
    case   9
        
        % ALG. DEG.:    9
        % PTS CARD.:   19
        % NEG. W.  :    0
        % OUT PTS. :    0
        % M.E.INF. : 3.1e-16
        
        xyw_bar=[
            3.33333333333335257719909350271337e-01 3.33333333333331482961625624739099e-01 3.33333333333333259318465024989564e-01 9.71357962827988224985276133338630e-02
            4.89682519198736176946340492577292e-01 4.89682519198739063526204517984297e-01 2.06349616025247595274549894384108e-02 3.13347002271390229211078803928103e-02
            4.89682519198737009613608961444697e-01 2.06349616025245374828500644071028e-02 4.89682519198738397392389742890373e-01 3.13347002271387731209273397325887e-02
            2.06349616025248878969922117221358e-02 4.89682519198738230858936049116892e-01 4.89682519198736898591306498929043e-01 3.13347002271393143546518444964022e-02
            4.37089591492937745709213004374760e-01 4.37089591492935525263163754061679e-01 1.25820817014126673516472010305733e-01 7.78275410047742782770896496913338e-02
            4.37089591492937745709213004374760e-01 1.25820817014125674315749847664847e-01 4.37089591492936524463885916702566e-01 7.78275410047739035768188387010014e-02
            1.25820817014127700472769788575533e-01 4.37089591492935580774314985319506e-01 4.37089591492936690997339610476047e-01 7.78275410047746946107238841250364e-02
            1.88203535619033524017851277676527e-01 1.88203535619031969705616802457371e-01 6.23592928761934617298834382381756e-01 7.96477389272102626049942841746088e-02
            1.88203535619033940351485512110230e-01 6.23592928761932951964297444646945e-01 1.88203535619033135439792658871738e-01 7.96477389272104707718114013914601e-02
            6.23592928761936171611068857600912e-01 1.88203535619031553371982568023668e-01 1.88203535619032275016948574375419e-01 7.96477389272100266826015513288439e-02
            4.47295133944526565605848134055123e-02 4.47295133944527745217811798283947e-02 9.10540973211094617489891334116692e-01 2.55776756586980312524470804191878e-02
            4.47295133944526010494335821476852e-02 9.10540973211094728512193796632346e-01 4.47295133944526357439031016838271e-02 2.55776756586979271690385218107622e-02
            9.10540973211094395445286409085384e-01 4.47295133944528369718263149934501e-02 4.47295133944527675828872759211663e-02 2.55776756586981388053025909812277e-02
            3.68384120547366744613526634566369e-02 7.41198598784498008384957756788936e-01 2.21962989160765289398113964125514e-01 4.32835393772897300546098620088742e-02
            3.68384120547364316000660267036437e-02 2.21962989160765927776353123590525e-01 7.41198598784497564295747906726319e-01 4.32835393772896051545195916787634e-02
            7.41198598784497897362655294273281e-01 3.68384120547358973052354258470587e-02 2.21962989160766205332109279879660e-01 4.32835393772890292263255673788080e-02
            7.41198598784497564295747906726319e-01 2.21962989160766371865562973653141e-01 3.68384120547360638386891196205397e-02 4.32835393772892443320365885028878e-02
            2.21962989160765400420416426641168e-01 3.68384120547361401665220626000519e-02 7.41198598784498452474167606851552e-01 4.32835393772891610653097416161472e-02
            2.21962989160765011842357807836379e-01 7.41198598784498452474167606851552e-01 3.68384120547364801723233540542424e-02 4.32835393772895010711110330703377e-02
            ];
        
    case  11
        
        % ALG. DEG.:   11
        % PTS CARD.:   28
        % NEG. W.  :    0
        % OUT PTS. :    0
        % M.E.INF. : 4.2e-16
        
        xyw_bar=[
            3.33333333333333425851918718763045e-01 3.33333333333333370340767487505218e-01 3.33333333333333148296162562473910e-01 8.79773011622321937652557721776247e-02
            2.59891409282874893960091355893383e-02 2.59891409282874824571152316821099e-02 9.48021718143424951819042689749040e-01 8.74431155373607202352381762011646e-03
            2.59891409282874408237518082387396e-02 9.48021718143425062841345152264694e-01 2.59891409282875240904786551254801e-02 8.74431155373604947211862992162423e-03
            9.48021718143425173863647614780348e-01 2.59891409282874408237518082387396e-02 2.59891409282873853126005769809126e-02 8.74431155373604773739515394481714e-03
            9.42875026479226829856372660287889e-02 9.42875026479227384967884972866159e-02 8.11424994704154634028725467942422e-01 3.80815719939349636713465940829337e-02
            9.42875026479225580855469956986781e-02 8.11424994704154745051027930458076e-01 9.42875026479227384967884972866159e-02 3.80815719939349636713465940829337e-02
            8.11424994704154634028725467942422e-01 9.42875026479225997189104191420483e-02 9.42875026479227662523641129155294e-02 3.80815719939349428546648823612486e-02
            4.94636775017213814464867027709261e-01 4.94636775017213869976018258967088e-01 1.07264499655723155591147133236518e-02 1.88554480561312423625430767515354e-02
            4.94636775017213758953715796451434e-01 1.07264499655723415799668529757582e-02 4.94636775017213869976018258967088e-01 1.88554480561312527708839326123780e-02
            1.07264499655723502535842328597937e-02 4.94636775017213814464867027709261e-01 4.94636775017213814464867027709261e-01 1.88554480561312597097778365196064e-02
            2.07343382614511434480775164956867e-01 2.07343382614511434480775164956867e-01 5.85313234770977075527298438828439e-01 7.21596975447395400093952844144951e-02
            2.07343382614511378969623933699040e-01 5.85313234770977297571903363859747e-01 2.07343382614511351214048318070127e-01 7.21596975447395122538196687855816e-02
            5.85313234770977297571903363859747e-01 2.07343382614511351214048318070127e-01 2.07343382614511351214048318070127e-01 7.21596975447395122538196687855816e-02
            4.38907805700492092970677049379447e-01 4.38907805700492148481828280637274e-01 1.22184388599015814058645901241107e-01 6.93291387055358782065539458017156e-02
            4.38907805700492203992979511895101e-01 1.22184388599015814058645901241107e-01 4.38907805700491926437223355605965e-01 6.93291387055358782065539458017156e-02
            1.22184388599015911203160555942304e-01 4.38907805700492148481828280637274e-01 4.38907805700491926437223355605965e-01 6.93291387055359059621295614306291e-02
            7.85647334377654633569776739674958e-17 8.58870281282636427455656757956604e-01 1.41129718717363461522040779527742e-01 7.36238378330058337167818294233257e-03
            1.57669518850046390999681416911743e-16 1.41129718717363572544343242043396e-01 8.58870281282636316433354295440949e-01 7.36238378330061546406248851326382e-03
            8.58870281282636427455656757956604e-01 1.42062476402819843287610091505200e-16 1.41129718717363433766465163898829e-01 7.36238378330060505572163265242125e-03
            8.58870281282636427455656757956604e-01 1.41129718717363350499738317012088e-01 2.22044604925031308084726333618164e-16 7.36238378330062066823291644368510e-03
            1.41129718717363572544343242043396e-01 2.02916231519083255217512748366332e-16 8.58870281282636205411051832925295e-01 7.36238378330063020921203431612412e-03
            1.41129718717363378255313932641002e-01 8.58870281282636427455656757956604e-01 2.22044604925031308084726333618164e-16 7.36238378330061199461553655964963e-03
            4.48416775891304561496575331602799e-02 6.77937654882590168270439789921511e-01 2.77220667528279340885433157382067e-01 4.10563154292885590379569293872919e-02
            4.48416775891304908441270526964217e-02 2.77220667528279451907735619897721e-01 6.77937654882590057248137327405857e-01 4.10563154292885729157447372017486e-02
            6.77937654882590168270439789921511e-01 4.48416775891304561496575331602799e-02 2.77220667528279396396584388639894e-01 4.10563154292885590379569293872919e-02
            6.77937654882590168270439789921511e-01 2.77220667528279340885433157382067e-01 4.48416775891304908441270526964217e-02 4.10563154292885451601691215728351e-02
            2.77220667528279340885433157382067e-01 4.48416775891305255385965722325636e-02 6.77937654882590168270439789921511e-01 4.10563154292885798546386411089770e-02
            2.77220667528279285374281926124240e-01 6.77937654882590279292742252437165e-01 4.48416775891303798218245901807677e-02 4.10563154292885590379569293872919e-02
            ];
        
end





