
function xyzw=cub_sphere_QMC(N,center,radius,pointset_type)

%-----------------------------------------------------------------------
% OBJECT
%-----------------------------------------------------------------------
% This routine determines a QMC cubature rule on the (outer) surface of the
% sphere defined by "center" 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).
%-----------------------------------------------------------------------
% INPUT
%-----------------------------------------------------------------------
% center: 1 x 3 matrix describing the centers in cartesian coordinates;
% radius: scalar, corresponding radius of the sphere;
% N: number of points in the 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 coor. and xyzw(:,4) are the weights of the QMC rule.
%-----------------------------------------------------------------------
% 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: November 6, 2024
%-----------------------------------------------------------------------


% Troubleshooting
if nargin < 1, N=1000; end % Points on the sphere 
if nargin < 2, center=[0 0 0]; end % Points on the sphere 
if nargin < 3, radius=1; end % Points on the sphere 
if nargin < 4, pointset_type = 'H'; end % QMC type

% 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

% 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 = [sqrt(1-X(:,1).^2).*cos(X(:,2)), ...
    sqrt(1-X(:,1).^2).*sin(X(:,2)), ...
    X(:,1)];

W=4*pi*radius^2*ones(size(S,1),1)/N;
xyzw=[S W];


