
function [IH,IL,dbox,flag,iters]=cub_adaptive(domain_structure,f,tol)

%--------------------------------------------------------------------------
% Object:
% Define algebraic cubature rule of algebraic degree of precision "deg" on
% domain defined by "domain_structure".
%--------------------------------------------------------------------------
% Input:
% domain_structure: structure defining the domain.
% f: function to integrate;
% tol: cubature tolerance
%--------------------------------------------------------------------------
% XW: cubature rule of degree "deg".
%--------------------------------------------------------------------------
% Example of "domain_structure":
%
% Polygon:
%      domain_struct.domain='polygon';
%      domain_struct.vertices=[0.1 0; 0.7 0.2; 1 0.5; 0.75 0.85];
% It is a polygon with vertices described by rows of the value of
% "domain_struct.vertices". The polygon is described counterclockwise by
% vertices.
%--------------------------------------------------------------------------

% ......................... troubleshooting ...............................

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

if ~isempty(domain_structure)

    % ................... check input ...................

    if ~isa(domain_structure,'struct')
        error('MATLAB: "define_cub_rule"',...
            getString(message(...
            'Variable "domain_structure" not defined')));
    end

    % ................... define domain ...................

    if isfield(domain_structure,'domain')
        domain = domain_structure.domain;
    else
        error('MATLAB: "define_cub_rule"',...
            getString(message(...
            'Variable "domain" not defined')));
    end

    % ..... domain parms check ......

    S(1)=strcmp(domain,'polygon');
    %     S(2)=strcmp(domain,'lune');
    %     S(3)=strcmp(domain,'circular-annular-sector');
    %     S(4)=strcmp(domain,'disk');
    %     S(5)=strcmp(domain,'asymmetric-circular-sector');
    %     S(6)=strcmp(domain,'asymmetric-annulus');
    %     S(7)=strcmp(domain,'vertical-circular-zone');
    %     S(8)=strcmp(domain,'horizontal-circular-zone');
    %     S(9)=strcmp(domain,'circular-segment');
    %     S(10)=strcmp(domain,'symmetric-lens');
    %     S(11)=strcmp(domain,'butterfly');
    %     S(12)=strcmp(domain,'candy');
    %     S(13)=strcmp(domain,'convex-level-curve');
    %     S(14)=strcmp(domain,'circle-arc');
    %     S(15)=strcmp(domain,'sphere');
    %     S(16)=strcmp(domain,'spherical-rectangle');
    %     S(17)=strcmp(domain,'pyramid');
    %     S(18)=strcmp(domain,'3D_rect');
    %     S(19)=strcmp(domain,'rotation-domain');
    %     S(20)=strcmp(domain,'polygonal-boundary');
    %     S(21)=strcmp(domain,'pluriinterval');
    %     S(22)=strcmp(domain,'sector');
    %     S(23)=strcmp(domain,'pyramid-boundary');
    %     S(24)=strcmp(domain,'asymmetric-circular-sector-boundary');
    %     S(25)=strcmp(domain,'asymmetric-annulus-boundary');
    %     S(26)=strcmp(domain,'vertical-circular-zone-boundary');
    %     S(27)=strcmp(domain,'horizontal-circular-zone-boundary');
    %     S(28)=strcmp(domain,'circular-segment-boundary');
    %     S(29)=strcmp(domain,'symmetric-lens-boundary');
    %     S(30)=strcmp(domain,'butterfly-boundary');
    %     S(31)=strcmp(domain,'candy-boundary');
    %     S(32)=strcmp(domain,'sector-boundary');
    %     S(33)=strcmp(domain,'rotation-surface');
    %     S(34)=strcmp(domain,'spherical-triangle');

end


% .......................... WAM definition ...............................
dbox=[];

switch domain

    case {'polygon','square','rectangle','triangle','unit-simplex',...
            'unit-square[0,1]x[0,1]'};

        vertices=domain_structure.vertices;
        [IH,IL,flag,iters]=cub_polygon_adaptive(vertices,f,tol);
        S=class(vertices);
        if strcmp(S,'polyshape')
            % ... XY is a polyshape object ...
            [xLimit, yLimit] = boundingbox(vertices);
            dbox=[xLimit; yLimit]';
        else
            % ... XY is a matrix ...
            x=vertices(:,1); mx=min(x); Mx=max(x);
            y=vertices(:,2); my=min(y); My=max(y);
            dbox=[mx Mx; my My]';
        end

    case 'spherical-triangle'
        sphtri_vertices=domain_structure.vertices;
        XW = cub_sphtri(deg,sphtri_vertices(1,:),sphtri_vertices(2,:),...
            sphtri_vertices(3,:));
        dbox=[];

    case 'spherical-polygon'
        sphpoly_vertices=domain_structure.vertices;
        [IL,Ierr,flag,IH,iters,tri_vertices,tri_conn_list]=...
            cub_sphpgon_adaptive(sphpoly_vertices,f,tol);
        dbox=[];


end






