
function [I,Ierr,flag,Ihigh,iters]=cub_sphrect_adaptive_2020(sphrect_intv,...
    f,tol)

%--------------------------------------------------------------------------
% Input:
% sphrect_intv: matrix 2 x 2;
%     Horizontally, "sphrect_intv" must denote the start point and the end 
%     point for each range.
%     Notice that 
%     * the first interval I1 must be contained in [0,pi], 
%     * the second interval I2 must be contained in [0,2*pi])
%     and that in cartesian coordinates
%
%           x_1=sin(theta_1)*cos(theta_2),
%           x_2=sin(theta_1)*sin(theta_2),
%           x_3=cos(theta_1),
%
%     with "theta1" in "I1" and "theta2" in "I2";
%
% f        : function to integrate over the spherical rectangle defined by
%            "sphrect_intv";
% tol      : absolute and relative cubature error tolerance; the stopping
%      criterion is of the form:
%
%              (abserr < tol) | abserr > tol*(abs(I)+eps)
%--------------------------------------------------------------------------
% Output:
% I      : approximation of integral of "f" over the spherical rectangle
%          defined by vertices
%--------------------------------------------------------------------------
% Subroutines:
% 1. generate_sphrect_sons.
%--------------------------------------------------------------------------

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

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

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

% rule degree of precision
ade=4;

% ........ Initial spherical 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 spherical 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 spherical 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={sphrect_intv};
[~,xyzw]=cub_sphrect(sphrect_intv,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.sphrect_intv : 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 
%                  spherical coordinates;
% * L2_data.integral : scalar representing the integral of "f" on the j-th
%           son sph. rectangle.
%..........................................................................

L2_data=generate_sphrect_sons(sphrect_intv,f,ade);
L2_data={L2_data};

% ... Evaluate absolute error ...
L2_data_temp=L2_data{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) > tol) & (sum(L1_errors) > tol*(abs(Itemp_high)+eps))
    
    N=length(L1_integrals);
    
    % too many rectangles: exit with errors
    if N > max_sphrect
        I=sum(L1_integrals);
        Ierr=sum(L1_errors);
        flag=1;
        if nargout > 3, Ihigh=Ihigh_computation(L2_data); end
        return;
    end
    
    [Ierr_max,kmax]=max(L1_errors);
    
    % Erase "kmax" sph. rectangle from LIST 1
    k_ok=setdiff(1:N,kmax);
    
    if length(k_ok) > 0
        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 length(k_ok) > 0
        % 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_sphrect_sons(L2_data_temp.sphrect_intv,...
            f,ade);
        
        % "L2_data_kmax_j_sons" is a cell with 4 structs data.
        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.sphrect_intv;
        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);
flag=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 spherical
%          rectangle.
%--------------------------------------------------------------------------

M=length(L2_data);
Ihigh=0;
for j=1:M
    L2_data_temp=L2_data{j};
    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_sphrect_sons(sphrect_intv,f,ade)

%--------------------------------------------------------------------------
% Input:
% sphrect_intv: matrix 2 x 2;
%     Horizontally, "sphrect_intv" must denote the start point and the end 
%     point for each range.
%     Notice that 
%     * the first interval I1 must be contained in [0,pi], 
%     * the second interval I2 must be contained in [0,2*pi])
%     and that in cartesian coordinates
%
%           x_1=sin(theta_1)*cos(theta_2),
%           x_2=sin(theta_1)*sin(theta_2),
%           x_3=cos(theta_1),
%
%     with "theta1" in "I1" and "theta2" in "I2";
%
% f  : function to integrate over the spherical rectangle defined by
%      vertices;
%
% ade: rule algebraic degree of precision
%--------------------------------------------------------------------------
% Output:
% L2_data :
%--------------------------------------------------------------------------

if nargin < 3, ade=4; end

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

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

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









function L2_dataL=make_L2_data(sphrect_intv,f,ade)

%--------------------------------------------------------------------------
% Input:
% sphrect_intv: matrix 2 x 2;
%     Horizontally, "sphrect_intv" must denote the start point and the end 
%     point for each range.
%     Notice that 
%     * the first interval I1 must be contained in [0,pi], 
%     * the second interval I2 must be contained in [0,2*pi])
%     and that in cartesian coordinates
%
%           x_1=sin(theta_1)*cos(theta_2),
%           x_2=sin(theta_1)*sin(theta_2),
%           x_3=cos(theta_1),
% f        : function to integrate over the spherical rectangle defined by
%            vertices.
% ade: rule algebraic degree of precision
%--------------------------------------------------------------------------
% Output:
% L2_dataL : determine from the spherical rectangle defined by vertices, a
%            struct data, with form:
%            * L2_dataL.sphrect_intv,
%            * L2_dataL.integral.
%--------------------------------------------------------------------------

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




