
function [I,Ierr,success,Ihigh,iters]=cub_torus_rect_adaptive(r,R,...
    omega,f,tol)

%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% This routine is a basic adaptive code, to compute quadrature over
% toroidal rectangles of continuous bivariate functions, up to a fixed 
% tolerance.
%
% As it is well-known, in specific and unfortunate situations, it may not
% provide correct results.
%
% The torus must not have self-intersections, that is equivalent to 
% "R > r".
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% r, R: minor and major radius of the torus (R < r),
%
% omega: matrix 2 x 2;
%     Horizontally, "omega" must denote the start point and the end 
%     point for each range.
%     Notice that 
%     * the first interval I1 must be contained in [-pi,pi], 
%     * the second interval I2 must be contained in [-pi,pi])
%     and that in cartesian coordinates
%
%           x(theta1,theta2)=(R+r*cos(theta1)).*cos(theta2);
%           y(theta1,theta2)=(R+r*cos(theta1)).*sin(theta2);
%           z(theta1,theta2)=r*sin(theta1);
%
%     with "theta1" in "I1" and "theta2" in "I2";
%
% f        : function to integrate over the toroidal rectangle defined by
%            "omega";
%
% tol   : vector of two components establishing the cubature error
%         tolerance
%
%            tol(1): 'AbsTol', absolute error tolerance
%            tol(2): 'RelTol', relative error tolerance
%
%         it attempts to satisfy |IH - I| <= max(AbsTol,RelTol*|IH|),
%         where "IH" is the numerical value computed by the algorithm;
%
%         if "tol" is scalar then "AbsTol=tol", "RelTol=tol".
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% I      : approximation of integral of "f" over the toroidal rectangle
%          defined by vertices;
%
% Ierr   : integration error;
%
% success: 0: right termination, 1: some issues;
%
% Ihigh  : fine approximation of the integral given by the algorithm;
%
% iters  : iterations of the algorithm.
%--------------------------------------------------------------------------
% Examples:
%--------------------------------------------------------------------------
% >>
% >> 
%--------------------------------------------------------------------------
% Subroutines:
%--------------------------------------------------------------------------
% 1. generate_torical_sons.
% Required routines:
%
% This program will use:
%
% 1. trigauss
%--------------------------------------------------------------------------
%% Copyright (C) 2026 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 <alvise@math.unipd.it>
%%          Marco Vianello   <marcov@math.unipd.it>
%%
%% Date       : January 6, 2026
%% Last Update: January 6, 2026
%--------------------------------------------------------------------------

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

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

% max number of sph. rectangles in which the domain is partitioned
max_sphrect=1000;

% rule degree of precision
ade=8;

% safety parameter (enforcing error estimate)
safe_parm=0.1;

% ........ Initial toroidal rectangle data stored in structure ............

%..........................................................................
% Main strategy:
%
% Below, we make 2 lists of sph. rectangle data:
% 1. LIST 1: all the rectangles provide contribution to the integration
%    result and for each rectangle an error estimate is available;
% 2. LIST 2: all the rectangles are "sons" rectangles of some rectangle in
%    LIST 1; integrals are computed in each of them but no error estimate
%    is available.
%
% We store the relevant data of each toroidal rectangle (LIST 1 set).
%
% * L1_rectangles: cell, whose k-th element are the ranges of the k-th
%                 sph. rectangle stored as a 2 x 2 matrix, whose j-th row
%                 are the range of the j-th toroidal coordinates;
% * L1_integrals: vector, whose k-component represent the approximation of
%                 the integral of "f" on the k-th sph. rectangle;
% * L1_errors   : vector, whose k-component represent the approximation of
%                 the error of integral on the k-th sph. rectangle.
%..........................................................................



% vertices
L1_rectangles={omega};
[~,xyzw]=cub_torus_rect(r,R,omega,ade);
Itemp_low=(xyzw(:,4))'*feval(f,xyzw(:,1),xyzw(:,2),xyzw(:,3));
L1_integrals=Itemp_low; % approximation of integral on sph.rect.(few evals)

%..........................................................................
% Note:
% We subdivide each sph.rectangle in the list LIST 1 in 4 sph.rectangles, 
% that we call "sons".
%
% Of each son rectangle we make a structure:
%
% * L2_data.omega : ranges of the k-th sph. rectangle stored as a 
%                 2 x 2 matrix, whose j-th row is the range of the j-th 
%                  toroidal coordinates;
% * L2_data.integral : scalar representing the integral of "f" on the j-th
%           son sph. rectangle.
%..........................................................................

L2_data=generate_torical_sons(r,R,omega,f,ade);
L2_data={L2_data};

% ... Evaluate absolute error ...
L2_data_temp=L2_data{1};
L2_data_integrals=zeros(4,1);
for j=1:4, L2_data_integrals(j)=L2_data_temp{j}.integral; end

% better approximation of integral on rectangle (more evals)
Itemp_high=sum(L2_data_integrals);

L1_errors=abs(Itemp_high-Itemp_low);

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

iters=1;

while sum(L1_errors) > safe_parm*max(tol(1),tol(2)*abs(Itemp_high))
    
    N=length(L1_integrals);
    
    % too many rectangles: exit with errors
    if N > max_sphrect
        I=sum(L1_integrals);
        Ierr=sum(L1_errors);
        success=1;
        if nargout > 3, Ihigh=Ihigh_computation(L2_data); end
        return;
    end
    
    [~,kmax]=max(L1_errors);
    
    % Erase "kmax" sph. rectangle from LIST 1
    k_ok=setdiff(1:N,kmax);
    
    if not(isempty(k_ok))
        L1_rectangles={L1_rectangles{k_ok}};
        L1_integrals=L1_integrals(k_ok);
        L1_errors=L1_errors(k_ok);
    else
        L1_rectangles={};
        L1_integrals=[];
        L1_errors=[];
    end
    
    % Move sub sph. rectangles relative to "k-max" from LIST 2 to LIST 1.
    
    L2_data_kmax=L2_data{kmax}; % cell with 4 data structs
    
    if not(isscalar(k_ok))
        % erase cell relative to sub rectangles of kmax
        L2_data={L2_data{k_ok}};
    else
        L2_data={};
    end
    % rectangle, from LIST 2
    
    %......................................................................
    % Generate sons of each son of L2_data_kmax and compute approximate
    % integration errors for each son of L2_data_kmax.
    % 1. Each son of L2_data_kmax is moved to LIST 1.
    % 2. The son of each son of L2_data_kmax is moved to LIST 2.
    %......................................................................
    
    for j=1:4
        L2_data_temp=L2_data_kmax{j}; % struct data
        L2_data_temp_sons=generate_torical_sons(r,R,L2_data_temp.omega,...
            f,ade);
        
        % "L2_data_kmax_j_sons" is a cell with 4 structs data.
        L2_data_temp_sons_integral=zeros(4,1);
        for jj=1:4
            L2_data_temp_sons_integral(jj)=L2_data_temp_sons{jj}.integral;
        end
        Itemp_low=L2_data_temp.integral;
        Itemp_high=sum(L2_data_temp_sons_integral);
        Itemp_err=abs(Itemp_high-Itemp_low);
        
        % update LIST 1 with j-th sub rectangle "generated" by kmax.
        L1_rectangles{end+1}=L2_data_temp.omega;
        L1_integrals(end+1)=L2_data_temp.integral;
        L1_errors(end+1)=Itemp_err;
        
        % update LIST 2 with sons of j-th sub rectangle "generated" by kmax.
        L2_data{end+1}=L2_data_temp_sons;
    end
    
    iters=iters+1;
    
end

% ....................... providing results ...............................

I=sum(L1_integrals);
Ierr=sum(L1_errors);
success=0;

if nargout > 3, Ihigh=Ihigh_computation(L2_data); end








function Ihigh=Ihigh_computation(L2_data)

%--------------------------------------------------------------------------
% Object:
% Approximation of the integral in view of approximations in LIST 2.
%--------------------------------------------------------------------------
% Input:
% L2_data: struct defining LIST 2 sph. rectangles vertices, areas and
%          integral approximations.
%--------------------------------------------------------------------------
% Output:
% Ihigh:   approximation of integral of "f" over the initial toroidal
%          rectangle.
%--------------------------------------------------------------------------

M=length(L2_data);
Ihigh=0;
for j=1:M
    L2_data_temp=L2_data{j};
    L2_data_temp_sons_integral=zeros(4,1);
    for jj=1:4
        L2_data_temp_sons_integral(jj)=L2_data_temp{jj}.integral;
    end
    Ihigh=Ihigh+sum(L2_data_temp_sons_integral);
end











function L2_data=generate_torical_sons(r,R,omega,f,ade)

%--------------------------------------------------------------------------
% Input:
% omega: matrix 2 x 2;
%     Horizontally, "omega" must denote the start point and the end 
%     point for each range.
%     Notice that each point will be represented as
%
%           x(theta1,theta2)=(R+r*cos(theta1)).*cos(theta2);
%           y(theta1,theta2)=(R+r*cos(theta1)).*sin(theta2);
%           z(theta1,theta2)=r*sin(theta1);
%
%     where 
%
%       a) "theta1" belongs to [omega(1,1),omega(1,2)] that is 
%          contained in [-pi,pi].
%       b) "theta2" belongs to [omega(2,1),omega(2,2)] that is 
%          contained in [-pi,pi].
%
% f  : function to integrate over the toroidal rectangle defined by
%      vertices;
%
% ade: rule algebraic degree of precision
%--------------------------------------------------------------------------
% Output:
% L2_data :
%--------------------------------------------------------------------------

if nargin < 3, ade=4; end

xmin=omega(1,1); xmax=omega(1,2); xmid=(xmin+xmax)/2;
ymin=omega(2,1); ymax=omega(2,2); ymid=(ymin+ymax)/2;

% ........................ rectangle data .................................

L2_data{1}=make_L2_data(r,R,[xmin xmid; ymin ymid],f,ade); 
L2_data{2}=make_L2_data(r,R,[xmid xmax; ymin ymid],f,ade);
L2_data{3}=make_L2_data(r,R,[xmin xmid; ymid ymax],f,ade);
L2_data{4}=make_L2_data(r,R,[xmid xmax; ymid ymax],f,ade);









function L2_dataL=make_L2_data(r,R,omega,f,ade)

%--------------------------------------------------------------------------
% Input:
% omega: matrix 2 x 2;
%     Horizontally, "omega" must denote the start point and the end 
%     point for each range.
%     Notice that each point will be represented as
%
%           x(theta1,theta2)=(R+r*cos(theta1)).*cos(theta2);
%           y(theta1,theta2)=(R+r*cos(theta1)).*sin(theta2);
%           z(theta1,theta2)=r*sin(theta1);
%
%     where 
%
%       a) "theta1" belongs to [omega(1,1),omega(1,2)] that is 
%          contained in [-pi,pi].
%       b) "theta2" belongs to [omega(2,1),omega(2,2)] that is 
%          contained in [-pi,pi].
% f        : function to integrate over the toroidal rectangle defined by
%            vertices.
% ade: rule algebraic degree of precision
%--------------------------------------------------------------------------
% Output:
% L2_dataL : determine from the toroidal rectangle defined by vertices, a
%            struct data, with form:
%            * L2_dataL.omega,
%            * L2_dataL.integral.
%--------------------------------------------------------------------------

if nargin < 3, ade=4; end
L2_dataL.omega=omega;
[~,xyzw]=cub_torus_rect(r,R,omega,ade);
L2_dataL.integral=(xyzw(:,4))'*feval(f,xyzw(:,1),xyzw(:,2),xyzw(:,3));




