function [res,err_est]=set_adaptive(f,region,varargin)

%%%%%
%
% INPUT:
%   f:    Function to integrate
%   region: Domain of integration (toroidal rectangle)
%
%   Named Arguments:
%   tol:    Tolerance
%   maxit:  Maximum recursion level
%   err_method: Error estimation method
%           1: Higher order Gauss
%           2: Antigauss
%           3: Gauss-Kronrod
%   domain:
%           1: Torus
%           2: Sphere
%   coord: Coordinate system
%           0: Cartesian
%           1: Polar (spherical or toroidal)
%
% OUTPUT:
%   res:    Computed integral
%   err_est: Error estimate
%
%%%%%

%# define defaults at the beginning of the code so that you do not need to
%# scroll way down in case you want to change something or if the help is
%# incomplete
options = struct('tol',10^(-15),'maxit',5,'err_method',1,'domain',1,'coord',0,'plot',0);

%# read the acceptable names
optionNames = fieldnames(options);

%# count arguments
nArgs = length(varargin);
if round(nArgs/2)~=nArgs/2
   error('EXAMPLE needs propertyName/propertyValue pairs')
end

for pair = reshape(varargin,2,[]) %# pair is {propName;propValue}
   inpName = lower(pair{1}); %# make case insensitive

   if any(strmatch(inpName,optionNames))
      %# overwrite options. If you want you can test for the right class here
      %# Also, if you find out that there is an option you keep getting wrong,
      %# you can use "if strcmp(inpName,'problemOption'),testMore,end"-statements
      options.(inpName) = pair{2};
   else
      error('%s is not a recognized parameter name',inpName)
   end
end

err_method = options.('err_method');
maxit = options.('maxit');
tol = options.('tol');
domain = options.('domain');
coord = options.('coord');
plot = options.('plot');

format long e;
queue = [];
done = [];
counter = 0;

%tic;
% Computing the integral on the whole region first
[preres,preerr,evaluations] = order_adaptive(f,region,tol,err_method,domain,coord,1);
counter = counter + evaluations;
x = [region(1,:), region(2,:)];

queue = [queue; preres, 1, 0, x];
err_est = queue(end,2);

it=0;
while (size(queue,1) > 0 && err_est > tol)
   [m, index] = max(queue(:,2));
   rect = queue(index,:);
   flag_break = 0;

    % If we've already recursed maxit times,
    % return this result, don't split anymore
    %rect(3)
    if (rect(3) >= maxit)
            queue = [queue(1:index-1,:); queue(index+1:end,:)];
            done = [done; rect];
            flag_break = 1;
    else

        % Split the region in four subregions,
        % Split along both directions
        xs(:,1) = rect(4:end);
        xs(:,2) = rect(4:end);
        xs(:,3) = rect(4:end);
        xs(:,4) = rect(4:end);
        part_res = [];

        hdist = (xs(2,1)-xs(1,1))/2;
        vdist = (xs(4,1)-xs(3,1))/2;
        if (hdist < 10^(-16) || vdist < 10^(-16))
                queue = [queue(1:index-1,:); queue(index+1:end,:)];
                done = [done; rect];
                flag_break = 1;
        else

            xs(2,1) = xs(2,1) - hdist;
            xs(4,1) = xs(4,1) - vdist;
            xs(1,2) = xs(1,2) + hdist;
            xs(4,2) = xs(4,2) - vdist;
            xs(2,3) = xs(2,3) - hdist;
            xs(3,3) = xs(3,3) + vdist;
            xs(1,4) = xs(1,4) + hdist;
            xs(3,4) = xs(3,4) + vdist;

            % Halve the tolerance if it is not too small to do so
            if (tol > 10^(-16))
                tol = tol/4;
            end

            % Actually calculate the partial results on the two parts
            prequeue = [queue(1:index-1,:)];
            for i=1:4
                reg = xs(:,i)';
                [p_res,p_err,p_eval] = order_adaptive( ...
                    f,[reg(1:2);reg(3:end)],tol,err_method,domain,coord,1);
                counter = counter + p_eval;
                if p_err > tol
                    prequeue = [prequeue; p_res, p_err, rect(3)+1, reg];
                else
                    done = [done; p_res, p_err, rect(3)+1, reg];
                end
            end
            prequeue = [prequeue; queue(index+1:end,:)];
            queue = prequeue;
        end
    end

    res = 0;
    err_est = 0;
    if size(queue,1) >= 1
        res  = res + sum(queue(:,1));
        err_est  = err_est + sum(queue(:,2));
    end
    if size(done,1) >= 1
        res  = res + sum(done(:,1));
        err_est  = err_est + sum(done(:,2));
    end
    old_err = err_est;
    preres = res;

end
%toc;
printf('Valutazioni della funzione f: %d\n',counter);
err_est = old_err;

% Plot center points of subdomains (only in interactive mode)
if plot > 0
    scatter_plot = [];
    for j = 1:size(done,1)
        rect = done(j,4:end);
        mid_rect = [(rect(1)+rect(2))/2, (rect(3)+rect(4))/2];
        t = mid_rect(1); p = mid_rect(2);
        mid_rect_cart = [(2+cos(t))*cos(p),(2+cos(t))*sin(p),sin(t)];
        scatter_plot = [scatter_plot; mid_rect_cart];
    end
    for j = 1:size(queue,1)
        rect = queue(j,4:end);
        mid_rect = [(rect(1)+rect(2))/2, (rect(3)+rect(4))/2];
        t = mid_rect(1); p = mid_rect(2);
        mid_rect_cart = [(2+cos(t))*cos(p),(2+cos(t))*sin(p),sin(t)];
        scatter_plot = [scatter_plot; mid_rect_cart];
    end

    [u,v]=meshgrid(region(1,1):(region(1,2)-region(1,1))/20:region(1,2),region(2,1):(region(2,2)-region(2,1))/20:region(2,2));
    x=(2+cos(v)).*cos(u);
    y=(2+cos(v)).*sin(u);
    z=sin(v);
    hold on;
    colormap(0.5*ones(size(z,1),3));
    axis("off");
    mesh(x,y,z)
    plot3(scatter_plot(:,1),scatter_plot(:,2),scatter_plot(:,3),'r.');
    view(180,90);
    hold off;
end

end
