
function xyzw=cub_sphpoly_QMC(N,vertices,pointset_type)

%-----------------------------------------------------------------------
% OBJECT
%-----------------------------------------------------------------------
% This routine determines a QMC cubature rule on the spherical polygon with
% vertices described (in order) by V;
% 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
%-----------------------------------------------------------------------
% N: number of points in the spherical polygon;
% vertices: N x 3 matrix, vertices of the polygon in spherical coordinates.
% 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, vertices=[1 0 0; 0 1 0; 0 0 1]; end % Points on the sphere 
if nargin < 3, 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 spherical polygon.
[X,CC,area_cap,area_sphpgon] = pts_in_sphpol(N,pointset_type,vertices);

% computing weights
W=area_sphpgon/size(X,1)*ones(size(X,1),1);

% providing output
xyzw=[X W];











function [X,CC,area_cap,area_sphpgon] = pts_in_sphpol(N,pts_type,vertices)
%--------------------------------------------------------------------------
% OBJECT.
%--------------------------------------------------------------------------
% This function gives as output a set of N points inside a spherical
% polygon defined through its vertices.
%
% The procedure works as follows.
%
% a) It rotates the spherical polygon "S" to the North Pole by a suitable
% rotation of its centroid. Let "S1" br such rotated sph. polygon.
%
% b) Next, it considers the "minimal" spherical cap "C" containing that
% mapped spherical polygon "S1" and determines a set of N uniform points
% in "C".
%
% c) Vertices and points of "C" are mapped via the stereographic map to the
% plane and are considered only those remaining inside the planar polygon.
%
% d) Finally the points are rotated back to the original sph. poly. "S".
%--------------------------------------------------------------------------
% INPUT
%--------------------------------------------------------------------------
% N        : Cardinality of the points inside the minimal spherical cap.
% pts_type : 'H' Halton points, 'S' Sobol points, 'R' random points.
% vertices : Cartesian coordinates of the vertices of the sph. polygon.
%--------------------------------------------------------------------------
% OUTPUT
%--------------------------------------------------------------------------
% X        : cartesian coordinates of the points inside the sph. polygon.
% CC       : centroid of the region.
% area_cap,area_sphpgon: areas of a small cap containing the spherical
%            polygon and of the spherical polygon itself.
%--------------------------------------------------------------------------
% DATES
%--------------------------------------------------------------------------
% First version: March, 2023;
% Checked: April 5, 2023.
%--------------------------------------------------------------------------
% AUTHORS
%--------------------------------------------------------------------------
% G. Elefante, A. Sommariva, M. Vianello
%--------------------------------------------------------------------------
% PAPER
%--------------------------------------------------------------------------
% "Qsurf: compressed QMC integration on algebraic surfaces"
% G. Elefante, A. Sommariva and M. Vianello
%--------------------------------------------------------------------------
% RELEASE DATE
%--------------------------------------------------------------------------
% First version: March 2023
% Last update: March, 2023
%--------------------------------------------------------------------------
% .................... compute barycenter vertices ........................

R=norm(vertices(1,:)); vertices=vertices/R;
CC=mean(vertices); CC=CC/norm(CC);

% ................ rotation matrix centroid to north pole .................

[az,el,r] = cart2sph(CC(1),CC(2),CC(3));
phi=az; theta=pi/2-el;
cp=cos(phi); sp=sin(phi); ct=cos(theta); st=sin(theta);
R1=[ct 0 -st; 0 1 0; st 0 ct]; R2=[cp sp 0; -sp cp 0; 0 0 1];
rotmat=R1*R2; inv_rotmat=rotmat';

% ........................ rotate vertices to north pole...................

vertices_NP=(rotmat*vertices')';

% ..................... uniform points in the sphere cap...................

z_min = min(vertices_NP(:,3));

% ..................... decide suff. card in the cap ......................
% Computing area by rule over spherical polygons.
XW=cub_sphpgon(1,vertices);
area_sphpgon=sum(XW(:,4));
h=R-z_min; area_cap=2*pi*R*h;

card_cap=ceil(3*area_cap*N/area_sphpgon);


switch pts_type
    case 'H'
        p = haltonset(2);
        Xsq = net(p,card_cap);
    case 'S'
        p = sobolset(2);
        Xsq = net(p,card_cap);
    case 'R'
        Xsq = rand(card_cap,2);
end

Xsq = [z_min+(1-z_min)*Xsq(:,1),2*pi*Xsq(:,2)];

XcapNP = [sqrt(1-Xsq(:,1).^2).*cos(Xsq(:,2)), ...
    sqrt(1-Xsq(:,1).^2).*sin(Xsq(:,2)), Xsq(:,1)];

% ................... stereographic map from south pole ...................

% ....... vertices .......

XX_SP=vertices_NP(:,1); YY_SP=vertices_NP(:,2); ZZ_SP=vertices_NP(:,3);
rat=1./(1+ZZ_SP);

XX_SPm=rat.*XX_SP; YY_SPm=rat.*YY_SP; % vertices on the plane

% ....... points .......

rat2=1./(1+XcapNP(:,3));
X_SP=rat2.*XcapNP(:,1); % points on the plane
Y_SP=rat2.*XcapNP(:,2);


% ........... polyshape of projected polygon from South Pole ..............

PG = polyshape(XX_SPm,YY_SPm);

in = isinterior(PG,X_SP,Y_SP);

Xin = [X_SP(in),Y_SP(in)]; % points in the planar polygonal domain

% mapping points to the spherical polygon at the North Pole
XinS = [2*Xin(:,1)./(1+Xin(:,1).^2+Xin(:,2).^2),...
    2*Xin(:,2)./(1+Xin(:,1).^2+Xin(:,2).^2),...
    (1-Xin(:,1).^2-Xin(:,2).^2)./(1+Xin(:,1).^2+Xin(:,2).^2)];

% mapping points to the required spherical polygon
X = (inv_rotmat*XinS')';

% selecting the first N points, as required by user
X=X(1:N,:);




