
function XW=cub_sphpgon(n,vertices)

%--------------------------------------------------------------------------
% Input:
% n        : algebraic degree of exactness of the rule;
% vertices : it is a M x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex of the spherical
%            polygon (counterclockwise order).
%--------------------------------------------------------------------------
% Output:
% XW       : N x 4 matrix, where the XW(k,1:3) are the cartesian
%            coordinates of the k-th node, while XW(k,4) is the respective
%            weight.
%--------------------------------------------------------------------------
% Subroutines:
% 1. triangulate_sphpgon_tg (attached);
% 2. cub_sphtri (external) and its subroutines.
%--------------------------------------------------------------------------

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

if norm(vertices(1,:)-vertices(end,:))==0,vertices=vertices(1:end-1,:);end

% ........ Procedure start up, data stored in structures ............

%..........................................................................
% Main strategy:
% First we triangulate the spherical polygon, providing the vertices of the
% triangulation and its connectivity list (description of each spherical
% triangle via its vertices).
% Next we determine a rule for each spherical triangle.
%..........................................................................


%[tri_vertices,tri_conn_list]=triangulate_sphpgon(vertices);
[tri_vertices,tri_conn_list]=triangulate_sphpgon_tg(vertices);


XW=[];

for k=1:size(tri_conn_list,1)
    verticesL=tri_vertices((tri_conn_list(k,:))',:);
    P1=verticesL(1,:); P2=verticesL(2,:); P3=verticesL(3,:);
    XWL = cub_sphtri(n,P1,P2,P3); XW=[XW; XWL];
end









function [tri_vertices,tri_conn_list]=triangulate_sphpgon(vertices)

%--------------------------------------------------------------------------
% Object:
% In this routine we determine a triangulation of the spherical polygon
% with M vertices described in cartesian coordinates the M x 3 matrix
% "vertices".
% The points of the triangulation are stored in "tri_vertices", while the
% K x 3 connectivity matrix is described in "tri_conn_list".
% In particular the k-th row of "tri_conn_list" consists of the indices of
% vertices of the k-th spherical triangle, w.r.t. "tri_vertices".
%--------------------------------------------------------------------------

if size(vertices,1) == 3
    tri_vertices=vertices; tri_conn_list=[1 2 3]; return;
end

% .................... compute barycenter vertices ........................

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

% .................... rotation matrix to the 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 ................................

vertices_NP=(rotmat*vertices')';

% ....... map radially vertices to plane tangent to north pole ............

XX_NP=vertices_NP(:,1); YY_NP=vertices_NP(:,2); ZZ_NP=vertices_NP(:,3);

% ............ polyshape of projected polygon to North Pole ...............

PG = polyshape(XX_NP,YY_NP);

% ............ triangulation of projected polygon to North Pole ...........

tri = triangulation(PG);

tri_vert_NP=PG.Vertices;

Xtri=tri_vert_NP(:,1); Ytri=tri_vert_NP(:,2); Ztri=sqrt(1-Xtri.^2-Ytri.^2);

tri_vertices_NP=[Xtri Ytri Ztri];
tri_vertices=R*(inv_rotmat*tri_vertices_NP')';
tri_conn_list=tri.ConnectivityList;









function [tri_vertices,tri_conn_list]=triangulate_sphpgon_tg(vertices)

%--------------------------------------------------------------------------
% Object:
% In this routine we determine a triangulation of the spherical polygon
% with M vertices described in cartesian coordinates the M x 3 matrix
% "vertices".
% The points of the triangulation are stored in "tri_vertices", while the
% K x 3 connectivity matrix is described in "tri_conn_list".
% In particular the k-th row of "tri_conn_list" consists of the indices of
% vertices of the k-th spherical triangle, w.r.t. "tri_vertices".
%--------------------------------------------------------------------------

if size(vertices,1) == 3
    tri_vertices=vertices; tri_conn_list=[1 2 3]; return;
end

% .................... compute barycenter vertices ........................

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

% .................... rotation matrix to the 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 ................................

vertices_NP=(rotmat*vertices')';

% ....... map radially vertices to plane tangent to north pole ............

XX_NP=vertices_NP(:,1); YY_NP=vertices_NP(:,2); ZZ_NP=vertices_NP(:,3);
rat=1./ZZ_NP;

XX_NPm=rat.*XX_NP; YY_NPm=rat.*YY_NP; ZZ_NPm=ones(size(XX_NPm));

% ............ polyshape of projected polygon to North Pole ...............

PG = polyshape(XX_NPm,YY_NPm);

% ............ triangulation of projected polygon to North Pole ...........

tri = triangulation(PG);

tri_vertices_NPm=PG.Vertices;
tri_vertices_NPm=[tri_vertices_NPm ones(size(tri_vertices_NPm,1),1)];
rads=sqrt((tri_vertices_NPm(:,1)).^2+(tri_vertices_NPm(:,2)).^2+1);

% ......................... output data ...................................

tri_vertices_NP=tri_vertices_NPm./rads;
tri_vertices=R*(inv_rotmat*tri_vertices_NP')';
tri_conn_list=tri.ConnectivityList;
