
function demo_plot

theta1=[pi/3 pi/4 pi/3];
theta2=[0 pi/3 pi/4];

vertices=[
    sin(theta1(1))*cos(theta2(1)) sin(theta1(1))*sin(theta2(1)) cos(theta1(1))
    sin(theta1(2))*cos(theta2(2)) sin(theta1(2))*sin(theta2(2)) cos(theta1(2))
    sin(theta1(3))*cos(theta2(3)) sin(theta1(3))*sin(theta2(3)) cos(theta1(3))
    ];

[tri_vertices,tri_conn_list]=triangulate_sphpgon(vertices);
plot_s2('spherical-triangle',tri_vertices,[],[],'',1,tri_vertices,...
    tri_conn_list);
% hold on:

vertices_data=generate_sphtri_sons(vertices);
for k=1:4
    vertices_L=vertices_data{k};
    [tri_vertices,tri_conn_list]=triangulate_sphpgon(vertices_L);
    plot_s2('spherical-triangle',tri_vertices,[],[],'',1,tri_vertices,...
        tri_conn_list);
end

hold off;
pause


function vertices_data=generate_sphtri_sons(vertices)

%--------------------------------------------------------------------------
% Input:
% vertices : it is a 3 x 3 matrix, where the k-th row represents the
%            cartesian coordinates of the k-th vertex (counterclockwise);
% f        : function to integrate over the spherical triangle defined by
%            vertices;
%--------------------------------------------------------------------------
% Output:
% vertices_data
%--------------------------------------------------------------------------

OA=(vertices(1,:)); OB=(vertices(2,:)); OC=(vertices(3,:));
R=norm(OA);

% ........................ compute midpoints ..............................

OAB_mid=(OA+OB)/2; OAB_mid=R*OAB_mid/norm(OAB_mid);
OAC_mid=(OA+OC)/2; OAC_mid=R*OAC_mid/norm(OAC_mid);
OBC_mid=(OB+OC)/2; OBC_mid=R*OBC_mid/norm(OBC_mid);

% ........................ triangle data ..................................

vertices=[OA; OAB_mid; OAC_mid];

vertices_data{1}=[OA; OAB_mid; OAC_mid];
vertices_data{2}=[OAB_mid; OB; OBC_mid];
vertices_data{3}=[OBC_mid; OC; OAC_mid];
vertices_data{4}=[OAB_mid; OBC_mid; OAC_mid];






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".
%--------------------------------------------------------------------------


% .................... 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;


