function [res,err,evaluations]=order_adaptive(f,x,tol,err_method,domain,coord,evaluations)

%%%%%
%
% INPUT:
%   f:      Function to integrate
%   x:      Region 
%   tol:    Tolerance
%   err_method: Error estimation method
%           1: Higher order Gauss
%           2: Antigauss
%           3: Gauss-Kronrod
%   domain:
%           1: Torus
%           2: Sphere
%   coord: Coordinate system
%           0: Cartesian (default)
%           1: Polar (spherical or toroidal)
%   evaluations: 
%           1: Calculate how many integrand evaluations 
%               have been performed
%           0: Don't do it (default)
%
% OUTPUT:
%   res:    Computed integral
%   err:    Error estimate
% 
%%%%%

if nargin < 7
    evaluations = 0;
end
if nargin < 6
    coord = 0;
end
if nargin < 5
    domain = 1;
end
if nargin < 4
    err_method = 1;
end
if nargin < 3
    tol = 10^(-16)
end

n=2;
res(1)=0;

cubature = cellfun(@str2func,{'cubature_2_torus' 'cubature_2_subsphere_fast'},'UniformOutput',false);
% Maximum number of iterations
maxexp=6;

% If we can't reach the required tolerance,
% try doubling the exactness degree of the formula
for i=2:maxexp
    n=2^i;

    % Computing nodes and weights of the formula
    [tw,xw]=cubature{domain}(x,n,0);

    % Computing the integral
    %fx=feval(f,xw(:,1),xw(:,2),xw(:,3));
    if coord > 0
        fx=feval(f,tw(:,1),tw(:,2));
        w=tw(:,3);
    else
        fx=feval(f,xw(:,1),xw(:,2),xw(:,3));
        w=xw(:,4);
    end
    if evaluations > 0
        evaluations = evaluations + 1;
    end

    %w=xw(:,4);
    res(i)=w'*fx;

    % Estimate the error using the specified method
    % Shifting the result indices so the routine returns
    % Kronrod or antigauss result in the end
    if err_method > 1
        [tw,xw]=cubature{domain}(x,n,0,err_method);
        fx=feval(f,xw(:,1),xw(:,2),xw(:,3));
        w=xw(:,4);
        res(i-1) = res(i);
        res(i) = w'*fx;
    end
    err=abs(res(i)-res(i-1));
    if err_method == 2
        err = err/2;
    end

    % We're done if the error is below the threshold in input
    if err < tol
        res=res(end);
        evaluations = evaluations - 1;
        return;
    end
end

evaluations = evaluations - 1;
res=res(end);
