
function [am_sphtr,deg_am] = am_sphtri(n,Q,AM_factor,center)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% The routine computes an admissible mesh (AM) "xyz" of
% degree "n" on the spherical triangle "S" with vertices "P1", "P2", "P3" 
% each one stored as row of "Q". 
%
% The sphere has center in the point "center" and radius defined by the
% vertices.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% n: admissible mesh degree (on the spherical triangle "S");
% Q: spherical triangle vertices, ordered counterclockwise;
%       they must be a "3 x 3" matrix, where "Q(k,:)" is a
%       vector containing the cartesian coordinates of the k-th vertex.
% * AM_factor: choose "AM_factor >= 5" but avoid values larger than "10"
%    since AM is better but with too many points (default: 7).
% * center: row vector with cartesian coordinates of the center of the 
%     sphere (default: [0 0 0]).
%
% Items with "*" above, may not be assigned by the user. In case of doubts,
% set them as "[]".
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% am_sphtr : 3-column array of nodes cartesian coords of the admissible
%             mesh of degree "n" on spherical triangle "S".
% deg_am: degrees of AM in spherical triangles partitioning the domain S.
%--------------------------------------------------------------------------
% Subroutines:
%--------------------------------------------------------------------------
% 1. compute_AM (internal)
% 2. wamquadrangle (internal)
% 3. triangle_subdivisions (internal)
% 4. am_sphtri_basic (internal)
% 5. alpha_omega_starV (internal)
%--------------------------------------------------------------------------
% Reference paper:
%--------------------------------------------------------------------------
% "Near-optimal polynomial interpolation on spherical triangles",
%  Mediterranean Journal of Mathematics, 9, 2022, article 68.
%  A. Sommariva and M. Vianello,
%--------------------------------------------------------------------------
% Dates:
%--------------------------------------------------------------------------
% Written on 08/01/2021: A. Sommariva and M. Vianello.
% Modified on 31/12/25 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------
% License:
%--------------------------------------------------------------------------
% Copyright (C) 2021- Alvise Sommariva, Marco Vianello.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%
% Authors:
%
% Alvise Sommariva <alvise@math.unipd.it>
% Marco Vianello <marcov@math.unipd.it>
%
% Date: December 31, 2025
%--------------------------------------------------------------------------


% ............................ troubleshooting ............................

% AM degree: ceil(AM_factor*alpha_omega_star*n)
if nargin < 3, AM_factor=[]; end
if isempty(AM_factor), AM_factor=7; end

if nargin < 4, center=[0 0 0]; end
if isempty(center), center=[0 0 0]; end


% .......... examine spherical triangle with two partition methods ........

% working on unit sphere and then provide result on given sphere

Q0=[Q(:,1)-center(1) Q(:,2)-center(2) Q(:,3)-center(3)]; 
R0=norm(Q0(1,:)); Q=Q0/R0;

[am_sphtr1,deg_am1]=compute_AM(n,Q,1,AM_factor);

% If the domain is decomposed in triangles try a second method, otherwise
% directly compute the AM.
if deg_am1 > 1
    [am_sphtr2,deg_am2]=compute_AM(n,Q,2,AM_factor);
    
    cards=[size(am_sphtr1,1) size(am_sphtr2,1)];
    [~,imin]=min(cards);
    
    % decide best option
    if imin == 1
        am_sphtr=am_sphtr1; deg_am=deg_am1;
    else
        am_sphtr=am_sphtr2; deg_am=deg_am2;
    end
else
    am_sphtr=am_sphtr1; deg_am=deg_am1;
end


% ............. determine AM on the given spherical triangle ..............

am_sphtr=R0*am_sphtr;
am_sphtr=[center(1)+am_sphtr(:,1) center(2)+am_sphtr(:,2) ...
    center(3)+am_sphtr(:,3)]; 








function [am_sphtr,deg_am]=compute_AM(n,Q,subd_meth,AM_factor)

%--------------------------------------------------------------------------
% Input:
% n: admissible mesh degree (on the spherical triangle "S");
% Q: spherical triangle vertices, ordered counterclockwise;
%       they must be a "3 x 3" matrix, where "Q(k,:)" is a
%       vector containing the cartesian coordinates of the k-th vertex.
% AM_factor: choose "AM_factor >= 5" but avoid values larger than "10"
%    since AM is better but with too many points (default: 7).
%--------------------------------------------------------------------------
% Output:
% am_sphtr : 3-column array of nodes cartesian coords of the admissible
%             mesh of degree "n" on spherical triangle "S".
% deg_am  : am degree, depending on subtriangle of Q.
%--------------------------------------------------------------------------
% Reference paper:
% "Chebyshev-Dubiner norming grids and Fekete-like interpolation on
% spherical triangles",
% A. Sommariva and M. Vianello,
% Department of Mathematics, University of Padova, Italy.
%--------------------------------------------------------------------------
% Dates:
% Written on 12/02/21 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

% Using degree: ceil(AM_factor*alpha_omega_star*n);

if nargin < 4, AM_factor=7; end


[Qlist1,omega_starV1]=triangle_subdivisions(Q,subd_meth);
L=length(Qlist1);

am_sphtr=[];

deg_am=zeros(L,1);

for k=1:L
    QL=Qlist1{k};
    omega_starL=omega_starV1(k);
    alpha_omega_starL=alpha_omega_starV(omega_starL);
    deg_am(k)=ceil(AM_factor*alpha_omega_starL*n);
    am_sphtrL = am_sphtri_basic(deg_am(k),QL);
    am_sphtr=[am_sphtr; am_sphtrL];
end









function [am_sphtr,am_tri_bar] = am_sphtri_basic(n,Q,am_tri_bar)

%--------------------------------------------------------------------------
% Object:
% The routine computes an admissible mesh (AM) "xyz" of
% degree "n" on the spherical triangle "S" with vertices "P1", "P2", "P3",
% stored as k-th row of "Q".
%
% This routine is based on
% * computing a AM of suitable degree on the triangle [0 0; 1 0; 0 1];
% * mapping the AM to the the planar triangle T with vertices
%   "Q" using barycentric coordinates;
% * mapping the AM from "T" to the spherical triangle "S" with vertices
%   "Q".
%
% Note: the spherical triangle is obtained by considering any point P of a
% planar triangle T and mapping radially the generic segment OP to the
% surface of the sphere.
%--------------------------------------------------------------------------
% Input:
% n: admissible mesh degree (on the spherical triangle "S");
% Q: spherical triangle vertices, ordered counterclockwise;
%       they must be a "3 x 3" matrix, where "Q(k,:)" is a
%       vector containing the cartesian coordinates of the k-th vertex.
%--------------------------------------------------------------------------
% Output:
% am_sphtr : 3-column array of nodes cartesian coords of the admissible
%             mesh of degree "n" on spherical triangle "S".
% am_tri_bar : 3-column array of nodes cartesian coords of the admissible
%             mesh of degree "n" on planar triangle "T" (in barycentric
%             coordinates).
%--------------------------------------------------------------------------
% Routines called:
% 1. wamquadrangle (internal)
%--------------------------------------------------------------------------
% Data:
% First version: 08/02/2021 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

if nargin < 3, am_tri_bar=[]; end

% .......... AM on a reference triangle (baricentric coordinates) .........
if isempty(am_tri_bar) == 1
    deg=n;
    Qref=[0 0; 1 0; 0 1; 0 0];
    am_tri=wamquadrangle(deg,Qref);
    am_tri_bar=[am_tri 1-sum(am_tri,2)];
end

% .......... mapping the AM on the planar triangle T ......................
am_tri=am_tri_bar*Q; % k-th point defined by k-th row


% ....... mapping the AM on T radially on the spherical triangle S ........
X=am_tri(:,1); Y=am_tri(:,2); Z=am_tri(:,3);
norm_am_tri=sqrt(X.^2+Y.^2+Z.^2);
am_sphtr=[X./norm_am_tri Y./norm_am_tri Z./norm_am_tri];









function wam = wamquadrangle(deg,Q)

%--------------------------------------------------------------------------
% Object:
% The routine produces a polynomial Weakly Admissible Mesh for degree "deg"
% in a convex planar quadrangle, namely
%
%    ||p||_quadrangle leq C(deg)*||p||_WAM with C(deg)=O(log^2(deg))
%
% Note: a grid for triangles can be obtained by collapsing some vertices of
% Q.
%
% Note: AM at degree "m" can be obtained by a suitable degree "deg >= m".
%--------------------------------------------------------------------------
% Input:
% deg: polynomial degree
% Q  : 2-columns array of the convex quadrangle vertices (the k-th vertex
%      is described by the k-th row).
%--------------------------------------------------------------------------
% Output:
% wam: 2-columns array of mesh points coordinates
%--------------------------------------------------------------------------
% Data:
% First version: 01/09/2015 by A. Sommariva and M. Vianello.
% Modified: 15/02/2021 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

% FUNCTION BODY

% Chebyshev-Lobatto grid
j=(0:1:deg);
c=cos(j*pi/deg);

A=Q(1,:); B=Q(2,:); C=Q(3,:); D=Q(4,:);

[u,v]=meshgrid(c);

pts=0.25*((1-u(:)).*(1-v(:))*A+(1+u(:)).*(1-v(:))*B);
pts=pts+0.25*((1+u(:)).*(1+v(:))*C+(1-u(:)).*(1+v(:))*D);

% eliminating possible multiple nodes in degenerate cases
% (triangle, segment)
dig=12;
wam=singlepts(pts,dig);








function [spts,ind] = singlepts(pts,dig)

%--------------------------------------------------------------------------
% Object:
% Eliminates from an array of points "pts" those having "dig" digits
% in common with a subset of other points, keeping only one per subset.
%--------------------------------------------------------------------------

atol=100*eps;
m=10.^(floor(log10(abs(pts))-dig+1));
pts1=round(pts./m).*m;
pts1(abs(pts)<atol) = 0;
[spts,ind]=unique(pts1,'rows');









function [CC,r,omega_star]=circumcenter(Q)
%--------------------------------------------------------------------------
% Object:
% Center and radius of the triangle Q.
%--------------------------------------------------------------------------
% Input:
% Q: spherical triangle vertices, ordered counterclockwise;
%    it is a "3 x 3" matrix, where "Q(k,:)" is a vector containing the 
%    cartesian coordinates of the k-th vertex.
%--------------------------------------------------------------------------
% Output:
% CC: circumcenter of the triangle Q.
% r: circumradius of the triangle Q.
% omega_star: angle(CC,0,CC+r).
%--------------------------------------------------------------------------
% Reference paper:
% "Chebyshev-Dubiner norming grids and Fekete-like interpolation on
% spherical triangles",
% A. Sommariva and M. Vianello,
% Department of Mathematics, University of Padova, Italy.
%--------------------------------------------------------------------------
% Reference:
%    https://en.wikipedia.org/wiki/Circumscribed_circle
%--------------------------------------------------------------------------
% Dates:
% Written on 15/02/21 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

A=Q(1,:); B=Q(2,:); C=Q(3,:);
a=A-C;
b=B-C;

r=norm(a)*norm(b)*norm(a-b)/(2*norm(cross(a,b)));

t1=(norm(a))^2*b-(norm(b))^2*a;
t2=cross(a,b);
t3=norm(t2);
CC=C+cross(t1,t2)/(2*t3^2);

omega_star=asin(r);









function [Qlist_OK,omega_starV]=triangle_subdivisions(Q,subd_method)

%--------------------------------------------------------------------------
% Object:
% We subdivide, if necessary, the initial triangle "Q" in a sequence of
% triangles "Qlist_OK" with "omega_star" small and "Qlist_OK" triangulation
% of "Q".
%--------------------------------------------------------------------------
% Input:
% Q: spherical triangle vertices, ordered counterclockwise;
%       they must be a "3 x 3" matrix, where "Q(k,:)" is a
%       vector containing the cartesian coordinates of the k-th vertex.
% subd_method: 
%    1. the triangle is decomposed in 4 triangles via side midpoints
%    2. the triangle is decomposed in 2 triangles via side midpoint of the
%    longer side.
%--------------------------------------------------------------------------
% Output:
% Qlist_OK: cell with N components, where the k-th component contains 
%    the matrix describing the k-triangle, via its rows.
% omega_starV: vector with N components, where the k-th is the angle of the 
%    smaller cap containing the k-th triangle described by "Qlist_OK{k}".
%--------------------------------------------------------------------------
% Reference paper:
% "Chebyshev-Dubiner norming grids and Fekete-like interpolation on
% spherical triangles",
% A. Sommariva and M. Vianello,
% Department of Mathematics, University of Padova, Italy.
%--------------------------------------------------------------------------
% Dates:
% Written on 13/02/20 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

if nargin < 2, subd_method=2; end

max_wstar=pi/5; % alternative: pi/6

Qlist_OK={};
Qlist_KO={Q};

omega_starV=[];
while isempty(Qlist_KO) == 0
    
    Q0=Qlist_KO{1};
    [~,~,omega_star]=circumcenter(Q0);
    
    % update lists
    Qlist_KO=Qlist_KO(2:end);
    
    if omega_star <= max_wstar
        % add to Qlist_OK
        Qlist_OK{end+1}=Q0;
        omega_starV(end+1,1)=omega_star;
    else
        % subdivide and add to Qlist_KO
        A=Q0(1,:); B=Q0(2,:); C=Q0(3,:);
        
        
        switch subd_method
            
            case 1
                
                % subdivide in 4 triangles via side midpoints
                
                AB0=(A+B)/2; AB=AB0/norm(AB0);
                AC0=(A+C)/2; AC=AC0/norm(AC0);
                BC0=(B+C)/2; BC=BC0/norm(BC0);
                
                QL1=[A; AB; AC]; QL2=[AB; B; BC];
                QL3=[AB; BC; AC]; QL4=[AC; BC; C];
                Qlist_KO{end+1}=QL1; Qlist_KO{end+1}=QL2;
                Qlist_KO{end+1}=QL3; Qlist_KO{end+1}=QL4;
                
                
            case 2
                
                % subdivide in 2 via midpoint of larger side
                
                l(1)=norm(A-B);  l(2)=norm(B-C); l(3)=norm(C-A);
                
                [~,imax]=max(l);
                
                switch imax
                    case 1
                        D=(A+B)/2; D=D/norm(D);
                        QL1=[A; D; C]; QL2=[B; C; D];
                    case 2
                        D=(B+C)/2; D=D/norm(D);
                        QL1=[A; B; D]; QL2=[A; D; C];
                    case 3
                        D=(C+A)/2; D=D/norm(D);
                        QL1=[A; B; D]; QL2=[B; C; D];
                end
                Qlist_KO{end+1}=QL1; Qlist_KO{end+1}=QL2;
                
                
        end
    end
end









function alpha_omega_star=alpha_omega_starV(omega_star)

%--------------------------------------------------------------------------
% Input:
% omega_star: angle(CC,0,CC+r) where CC is the circumcenter of the minimal
% cap containing the spherical triangle and r is its radius.
%--------------------------------------------------------------------------
% Output:
% alpha_omega_star: factor that is necessary to determine the degree of the
% AM to be used.
%--------------------------------------------------------------------------
% Reference paper:
% "Chebyshev-Dubiner norming grids and Fekete-like interpolation on
% spherical triangles",
% A. Sommariva and M. Vianello,
% Department of Mathematics, University of Padova, Italy.
%--------------------------------------------------------------------------
% Dates:
% Written on 12/02/21 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------


if omega_star <= pi/12, alpha_omega_star=1.05319e+00; return; end
if omega_star <= pi/10, alpha_omega_star=1.08000e+00; return; end
if omega_star <= pi/9, alpha_omega_star=1.10228e+00; return; end
if omega_star <= pi/8, alpha_omega_star=1.14; return; end
if omega_star <= pi/6, alpha_omega_star=1.31; return; end
if omega_star <= pi/5, alpha_omega_star=1.60; return; end
if omega_star <= 0.9*pi/4, alpha_omega_star=2.18; return; end
if omega_star <= 0.999*pi/4, alpha_omega_star=20.9; return; end
warning('Large omega_star');
alpha_omega_star=6.60204e+01;
