
function [xyzw,X,bbox,A] = qmc_union_spheres_26(centers,radii,N,pointset_type)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% This routine determines a QMC cubature rule on the (outer) surface of the
% union of spheres defined by there "centers" and "radii".
%
% In particular it provides a QMC cubature rule, stored in "xyzw", where
%
% 1. the nodes correspond to xyzw(:,1:3)
% 2. the weights correspond to xyzw(:,1:4)
%
% Furthermore, the value "A(k)" is the area of the subregion of the k-th
% sphere, that belongs to such (outer) surface.
%--------------------------------------------------------------------------
% INPUT
%--------------------------------------------------------------------------
% centers: N x 3 matrix describing the centers in cartesian coordinates;
%
% radii  : N x 1 vector, in which radii(k) is the radius of the k-th ball;
%
% N      : number of points in each sphere;
%
% pointset_type: string determining the 3D low discrepancy pointset to be
%          used:
%             'S': sobolset, 'H': haltonset, 'R': random
%--------------------------------------------------------------------------
% OUTPUT
%--------------------------------------------------------------------------
% xyzw: N x 4 matrix, in which xyzw(:,1:3) are the nodes described in
%        cartesian coordinates, xyzw(:,4) are the weights of the QMC rule.
%
% X   : N x 3 pointset in the bounding box, containing the nodes;
%
% bbox: bounding box as [xmin xmax; ymin ymax; zmin zmax]';
%
% A   : A(k) is the area of the subregion of the k-th sphere, that belongs 
%       to such (outer) surface.
%--------------------------------------------------------------------------
% AUTHORS
%--------------------------------------------------------------------------
% G. Elefante, A. Sommariva, M. Vianello
%--------------------------------------------------------------------------
% PAPER
%--------------------------------------------------------------------------
% "Qbubble: a numerical code for compressed QMC volume and surface
% integration on union of spheres"
% G. Elefante, A. Sommariva and M. Vianello
%--------------------------------------------------------------------------
%  RELEASE DATE
%--------------------------------------------------------------------------
% First version: November 2022
% Last update  : January 24, 2026 (A. Sommariva)
%--------------------------------------------------------------------------


% ........................... Troubleshooting .............................

% Points on each sphere

if nargin < 3
    Nsph=length(radii);
    N = ceil(2*10^6/Nsph);
end

if nargin < 4, pointset_type = 'H'; end




% ........................... Main code below .............................

% 1) Minimal box containing the spheres
a = min(centers-radii,[],1); b = max(centers+radii,[],1);



% 2) Choosing QMC points.
switch pointset_type
    case 'S'
        p = sobolset(2); X = net(p,N); clear p;
    case 'H'
        p = haltonset(2); X = net(p,N); clear p;
    case 'R'
        X = rand(N,2);
end



% 3) Defining points on the sphere.

% a) mapping points into the rectangle [-1,1]x[0,2pi]
X = [2*X(:,1)-1, X(:,2)*2*pi];

% b) mapping the points on the unitary sphere
S(:,1)=sqrt(1-X(:,1).^2).*cos(X(:,2));
S(:,2)=sqrt(1-X(:,1).^2).*sin(X(:,2));
S(:,3)=X(:,1);



% 4) Defining points on the union of spheres

Nsph = size(centers,1); % Number of centers
PtsBound = zeros(N*Nsph,3);

% Ordering Points in the all the spheres to obtain a sequence in the union
for i = 1:Nsph
    PtsBound(i:Nsph:end) = radii(i)*S + centers(i,:);
end



% Testing if the points are on the boundary/outside of each sphere or
% inside any other.

tol = -1e-14;

XN=( PtsBound(:,1)-(centers(:,1))').^2;
YN=( PtsBound(:,2)-(centers(:,2))').^2;
ZN=( PtsBound(:,3)-(centers(:,3))').^2;

Y=sum((XN+YN+ZN)-(radii.^2)' > tol,2);


% 5. Selecting the points on the boundary of the union

idx = (Y==Nsph);


% 6. Computing the area of the surface as percentage of each sphere and
%    QMC weights.

Yrsh=reshape(idx,Nsph,N);
card_sph_in=sum(Yrsh,2);
spheres_area=4*pi*(radii.^2);

A=(card_sph_in/N).*spheres_area;
wL0=A./card_sph_in;
wMAT=bsxfun(@times,Yrsh,wL0);
wV0w=wMAT(:);
idok=find(wV0w > 0);

xyz=PtsBound(idok,:);
ww=wV0w(idok,:);


% 7. Output matrix and a box containing the union of spheres

xyzw=[xyz ww];

if nargout > 2
    bbox=[a(1) b(1); a(2) b(2); a(3) b(3)]';
end



