function [xyz,deg_wam_2D] = mesh_sphtri_03(n,sphtri_vertices,toll)

%--------------------------------------------------------------------------
% Object:
% The routine computes a "numerical" weakly admissible mesh (WAM) "xyz" of
% degree "n" on the spherical triangle with vertices "P1", "P2", "P3".
%--------------------------------------------------------------------------
% Input:
% n: weakly admissible;
% sphtri_vertices: spherical triangle vertices, ordered counterclockwise;
%       they must be a "3 x 3" matrix, where "sphtri_vertices(k,:)" is a
%       vector containing the cartesian coordinates of the k-th vertex.
%--------------------------------------------------------------------------
% Output:
% xyz : 3-column array of nodes cartesian coords of the pointsed
%--------------------------------------------------------------------------
% Routines called:
% 1. compute_m (attached)
% 2. gqcircsect (external)
%--------------------------------------------------------------------------
% Data:
% First version: 08/01/2021 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

% ............................ troubleshooting ............................
P1=sphtri_vertices(1,:); P1=P1'; % column vectors
P2=sphtri_vertices(2,:); P2=P2';
P3=sphtri_vertices(3,:); P3=P3';

if nargin < 3, toll = 0; end

RR=norm(P1);


% ------------------------------- main code -------------------------------

% ... cicumcenter ...
% ... computing circumcenter ...
A=P1; B=P2; C=P3;
a=A-C; a_norm=norm(a);
b=B-C; b_norm=norm(b);
ab_cross=cross(a,b); ab_cross_norm=norm(ab_cross);
CC=cross((a_norm^2*b-b_norm^2*a),ab_cross)/(2*ab_cross_norm.^2)+C;

% .................... 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];
R=R1*R2; invR=R';

% ............... vertices of the triangle at the North Pole ..............

% A, B, C are column vectors
A=R*P1; B=R*P2; C=R*P3; vert=[A B C]; % A,B,C are column vectors

% ...... computing "m" needed in determining the degree of precision ......

% Using normalization factor RR on the projections A,B,C.
% In other words, we make up a WAM on the unit sphere (by
% contraction w.r.t. the unit sphere and map back the so obtained rule to
% the starting sphere with radius R).

% ... mapping vertices to tangent plane ...
dd(1)=1/A(3); dd(2)=1/B(3); dd(3)=1/C(3);
AA=dd(1)*A; BB=dd(2)*B; CC=dd(3)*C; % column vectors

vert0=[AA BB CC];

% ... wam on 2D planar triangle ....
t1=max(dd); m=compute_m(t1,toll);
deg_wam_2D=n*2*m;
wam_pts=wam_polygon(deg_wam_2D,(vert0(1:2,:))');

% ... wam on plane triangle on tangent plane North Pole ...
wam_ptsT=[wam_pts ones(size(wam_pts,1),1)];
rad_wam_ptsT=sqrt(wam_ptsT(:,1).^2+wam_ptsT(:,2).^2+wam_ptsT(:,3).^2);

% ... wam on sph. triangle on North Pole ...
wam_pts0=[wam_ptsT(:,1)./rad_wam_ptsT wam_ptsT(:,2)./rad_wam_ptsT ...
    wam_ptsT(:,3)./rad_wam_ptsT];

% ... wam on original sph. triangle ...
xyz=invR*wam_pts0'; xyz=xyz';

% Exporting results to the sphere with radius "R".
X=RR*xyz(:,1); Y=RR*xyz(:,2); Z=RR*xyz(:,3);
xyz=[X Y Z];









function wam_pts=wam_polygon(deg,polygon_vertices)

%--------------------------------------------------------------------------
% Object:
%
% The following routine computes a WAM over a polygon, of degree "deg".
%--------------------------------------------------------------------------
% Input:
%
% deg: WAM degree.
% polygon_vertices: polygon vertices, ordered counterclockwise; they must
%       be a "M x 2" matrix, where "polygon_vertices(k,:)" is a vector
%       containing the coordinates of the k-th vertex.
%--------------------------------------------------------------------------
% Output:
%
% wam_pts: wam points at degree "deg", it is a matrix "N x 2", where
%      "wam_pts(k,:)" is a vector containing the coordinates of the k-th
%      point of the wam.
%--------------------------------------------------------------------------
% Note:
% This routine requires the Matlab built-in "polyshape" environment.
%--------------------------------------------------------------------------
% Dates:
% Written on November 17, 2009, by F. Basaglia, A. Sommariva, M. Vianello.
%
% Modified on November 17, 2020, by A. Sommariva.
%--------------------------------------------------------------------------

% Determine polyshape object

xv=polygon_vertices(:,1); yv=polygon_vertices(:,2);
pgon = polyshape(xv,yv);

% Determine triangulation of polyshape object
tri = triangulation(pgon);
TCL=tri.ConnectivityList; N=size(TCL,1);
X = tri.Points(:,1); Y = tri.Points(:,2);


% Wam points in barycentric coordinates
[wam_pts,wam_pts_bar]=ref_wam_pts(deg);

% Wam points of the polygon
wam_pts=[];
for k=1:N
    
    % generating k-triangle of polygonal region triangulation
    TTloc=TCL(k,:);
    Xloc=X(TTloc); Yloc=Y(TTloc);
    
    % mapping reference points to triangle
    wam_pts_loc=wam_pts_bar*[Xloc Yloc];
    wam_pts=[wam_pts; wam_pts_loc];
end









function [wam_pts,wam_pts_bar]=ref_wam_pts(deg)

%--------------------------------------------------------------------------
% Object:
%
% This function determines a wam "wam" of degree "deg" on the reference
% simplex whose vertices are (0,0),(1,0),(0,1).
%--------------------------------------------------------------------------
% Input:
%
% M: WAM degree.
%--------------------------------------------------------------------------
% Output:
%
% wam_pts: points on the reference triangle (0,0),(1,0),(0,1); it is a
%      matrix "N x 2", where "wam_pts(k,:)" is a vector containing the
%      cartesian coordinates of the k-th  point of the wam.
%
% wam_pts_bar: points on the reference triangle (0,0),(1,0),(0,1); it is a
%      matrix "N x 3", where "wam_pts(k,:)" is a vector containing the
%      baricentric coordinates of the k-th  point of the wam.
%--------------------------------------------------------------------------

triangle=[0 1; 1 0; 0 0];

n1=2*deg;
n2=n1;
j1=(0:1:deg-1);
j2=(0:1:deg);
[rho,theta]=meshgrid(cos(j1*pi/n1),j2*pi/n2);

% Mapping to the simplex.

B=[rho(:).*cos(theta(:)) rho(:).*sin(theta(:))];

meshS=[B(:,1).^2 B(:,2).^2];
meshT(:,1)=triangle(3,1)*(1-meshS(:,1)-meshS(:,2))+...
    triangle(2,1)*meshS(:,2)+triangle(1,1)*meshS(:,1);
meshT(:,2)=triangle(3,2)*(1-meshS(:,1)-meshS(:,2))+...
    triangle(2,2)*meshS(:,2)+triangle(1,2)*meshS(:,1);

nodes_x=meshT(:,1); nodes_y=meshT(:,2);

nodes_x=[nodes_x; 0]; nodes_y=[nodes_y; 0];

wam_pts=[nodes_x nodes_y];
wam_pts_bar=[wam_pts 1-sum(wam_pts,2)];









function m=compute_m(t1,toll)

%--------------------------------------------------------------------------
% Object:
% It computes the "m" positive integer value depending from the sph. tri.
% with vertices P1, P2, P3, so that a numerical WAM rule over the sph.
% triangle of degree "n" can be obtained via a WAM of degree "m*n" on
% its projection on the xy-plane.
%--------------------------------------------------------------------------
% Input:
% t1: value to determine the interval [1,t1] of interest
%--------------------------------------------------------------------------
% Output:
% m: positive integer value depending from the sph.triangle with vertices
% P1, P2, P3, so that a numerical WAM rule over the sph.
% triangle of degree "n" can be obtained via a WAM of degree "m*n" on
% its projection on the xy-plane.
%--------------------------------------------------------------------------
% Data:
% First version: 08/01/2021 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

intvf=[1,t1]; F=@(x) 1/sqrt(x); f=chebfun(F,intvf);
m=max(1,(length(f)-1));
if toll > 0
    m0=m;
    err=0;
    while err <= toll
        [~,err]=minimax(f,m);
        m=m-1;
    end
        if m < m0, m=max(1,m+1); fprintf('\n \t REMEZ: IN: %3.0f OUT: %3.0f',m0,m); end
end