
function [I,Ierr,success,Ihigh,iters,tri_vertices,tri_conn_list]=...
    cub_sphpgon_adaptive(vertices,f,tol)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Adaptive cubature over spherical polygons contained in a cap whose
% polar angle is strictly inferior than π.
%
% The spherical polygon is on the sphere centered in the origin, whose 
% radius is defined by the "vertices" of the polygon.
%
% The integrand "f" must be given, as well as the tolerance "tol" of the
% absolute integration error.
%
% Note: this is an adaptive code and may fail in computing the desired
% integral with the required accuracy as all the adaptive codes.
%--------------------------------------------------------------------------
% Chebfun requirement:
%--------------------------------------------------------------------------
% This routine requires the installation of Chebfun.
% See "https://www.chebfun.org/download/" for details.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% 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);
% f        : function to integrate over the spherical polygon defined by
%            vertices;
% tol      : absolute cubature error tolerance.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% I      : approximation of integral of "f" over the spherical polygon
%           defined by vertices
% Ierr   : absolute error estimate
% success: 0: unsuccessful results, 1: successful results
% Ihigh  : integral on a finer triangulation
% tri_vertices: M x 3 matrix of vertices of the triangulation of the
%          spherical polygon into M triangles
% tri_conn_list: connectivity matrix (defines the sph.triangles in the 
%          triangulation) 
%--------------------------------------------------------------------------
% Subroutines:
%--------------------------------------------------------------------------
% 1. area_spherical_triangle;
% 2. center_spherical_triangle;
% 3. generate_sphtri_sons;
% 4. cub_sphtri.
%--------------------------------------------------------------------------
% Reference papers:
%--------------------------------------------------------------------------
% "Numerical cubature and hyperinterpolation over Spherical Polygons",
% Applied Mathematics and Computation,  Volume 495, 15 June 2025, 129335, 
%
% "Adaptive RBF cubature by scattered data on spherical polygons",
% R. Cavoretto, A. De Rossi, G. Elefante, A. Sommariva.
%--------------------------------------------------------------------------
%% 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 2 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, Marco Vianello.
%%
%% Date: JULY 07, 2021
%% Modified: DECEMBER 30, 2025
%--------------------------------------------------------------------------





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

if nargin < 3, tol=10^(-6); end
tri_meth=1; % Triangulation method.

% max number of triangles in which the domain is partitioned
max_triangles=1000;

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).
%
% Below, we make 2 lists of sph. triangle data:
% 1. LIST 1: all the triangles provide contribution to the integration
%    result and for each triangle an error estimate is available;
% 2. LIST 2: all the triangles are "sons" triangles of some triangle in
%    LIST 1; integrals are computed in each of them but no error estimate
%    is available.
%
% We store the relevant data of each spherical triangle (LIST 1 set).
%
% * L1_vertices : cell, whose k-th element are the vertices of the k-th
%                 sph. triangle stored as a 3 x 3 matrix, whose j-th row
%                 are the cartesian coordinates  of the j-th vertex of the
%                 k-th sph. triangle (ordered counterclockwise);
% * L1_areas    : vector, whose k-component represent the area of the k-th
%                 sph. triangle;
% * L1_integrals: vector, whose k-component represent the approximation of
%                 the integral of "f" on the k-th sph. triangle;
% * L1_errors   : vector, whose k-component represent the approximation of
%                 the error of integral on the k-th sph. triangle.
% 
% If the error estimate is not satisfied, the code find the "worst" 
% triangle, subdvides it, and updates the lists.
%..........................................................................


if tri_meth == 1
    [tri_vertices,tri_conn_list]=triangulate_sphpgon(vertices);
else
    [tri_vertices,tri_conn_list]=triangulate_sphpgon_tg(vertices);
end

L1_vertices={}; L1_integrals=[]; L1_errors=[]; L2_data={};


for k=1:size(tri_conn_list,1)
    vertices=tri_vertices((tri_conn_list(k,:))',:);
    [L1_vertices,L1_integrals,L1_errors,L2_data]=...
        start_up(vertices,f,L1_vertices,L1_integrals,L1_errors,L2_data);
end


% ........................ successive refinement ..........................

iters=1;

while  sum(L1_errors) > tol
    
    % fprintf('\n \t %6.0f',iters);
    
    N=length(L1_integrals);
    
    % too many triangles: exit with errors
    if N > max_triangles
        I=sum(L1_integrals); Ierr=sum(L1_errors); success=0;
        if nargout > 3, Ihigh=Ihigh_computation(L2_data); end
        return;
    end
    
    [~,kmax]=max(L1_errors);
    
    % Erase "kmax" sph. triangle from LIST 1
    k_ok=setdiff(1:N,kmax);
    
    if isempty(k_ok) 
        L1_vertices={L1_vertices{k_ok}}; L1_integrals=L1_integrals(k_ok);
        L1_errors=L1_errors(k_ok);
    else
        L1_vertices={}; L1_integrals=[]; L1_errors=[];
    end
    
    % Move sub sph. triangles relative to "k-max" from LIST 2 to LIST 1.
    
    L2_data_kmax=L2_data{kmax}; % cell with 4 data structs
    
    if isempty(k_ok)
        % erase cell relative to sub triangles of kmax
        L2_data={L2_data{k_ok}};
    else
        L2_data={};
    end
    % triangle, from LIST 2
    
    %......................................................................
    % Generate sons of each son of L2_data_kmax and compute approximate
    % integration errors for each son of L2_data_kmax.
    % 1. Each son of L2_data_kmax is moved to LIST 1.
    % 2. The son of each son of L2_data_kmax is moved to LIST 2.
    %......................................................................
    
    for j=1:4
        L2_data_temp=L2_data_kmax{j}; % struct data
        L2_data_temp_sons=generate_sphtri_sons(L2_data_temp.vertices,f);
        
        % "L2_data_kmax_j_sons" is a cell with 4 structs data.
        L2_data_temp_sons_integral=zeros(4,1);
        for jj=1:4
            L2_data_temp_sons_integral(jj)=L2_data_temp_sons{jj}.integral;
        end
        Itemp_low=L2_data_temp.integral;
        Itemp_high=sum(L2_data_temp_sons_integral);
        Itemp_err=abs(Itemp_high-Itemp_low);
        
        % update LIST 1 with j-th sub triangle "generated" by kmax.
        L1_vertices{end+1}=L2_data_temp.vertices;
        L1_integrals(end+1)=L2_data_temp.integral;
        L1_errors(end+1)=Itemp_err;
        
        % update LIST 2 with sons of j-th sub triangle "generated" by kmax.
        L2_data{end+1}=L2_data_temp_sons;
    end
    
    iters=iters+1;
    
end

% ....................... providing results ...............................

I=sum(L1_integrals); Ierr=sum(L1_errors); success=1;
if nargout > 3, Ihigh=Ihigh_computation(L2_data); 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".
%--------------------------------------------------------------------------


% .................... 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] = 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 [L1_vertices,L1_integrals,L1_errors,L2_data,Itemp_low,...
    Itemp_high]=start_up(vertices,f,L1_vertices,L1_integrals,L1_errors,...
    L2_data)

%--------------------------------------------------------------------------
% Object:
% Given a sph.triangle with vertices "vertices", a cell with previous
% vertices of previous triangles, a vector "L1_integrals" with integrals of
% "f" in the previous triangles and a vector "L1_errors" of the error
% estimates of the integrals of "f" in the previous triangles, it updates
% these cells and vectors with vertices of the new triangle, the integral
% of "f" in the new triangle, and its error estimate.
% Data about the sons of the new triangle is put in the cell "L2_data".
%--------------------------------------------------------------------------


% vertices
L1_vertices{end+1}=vertices;

% integrals
P1=(vertices(1,:))'; P2=(vertices(2,:))'; P3=(vertices(3,:))';

xyzw= cub_sphtri(4,P1,P2,P3);
Itemp_low=(xyzw(:,4))'*feval(f,xyzw(:,1),xyzw(:,2),xyzw(:,3));

% approximation of integral on triangle (few evals)
L1_integrals=[L1_integrals; Itemp_low];

%..........................................................................
% Note:
% We subdivide each sph. triangle in the list LIST 1 in 4 triangles, that
% we call "sons".
%
% Each initial planar triangle defining the generical sph. triangle, has
% 3 midpoints. In view of the midpoints, we get 4 triangles.
%
% Of each son triangle we make a structure:
%
% * L2_data.vertices : matrix 3 x 3, whose l-th row are the
%           cartesian coordinates of the l-th vertex of the j-th sub sph.
%           triangle;
% * L2_data.area     : scalar representing the area of the j-th son sph.
%           triangle;
% * L2_data.integral : scalar representing the integral of "f" on the j-th
%           son sph. triangle.
%..........................................................................

L2_data_temp=generate_sphtri_sons(vertices,f);
L2_data{end+1}=L2_data_temp;

% ... Evaluate absolute error ...
L2_data_integrals=zeros(4,1);
for j=1:4
    L2_data_integrals(j)=L2_data_temp{j}.integral; 
end

% better approximation of integral on triangle (more evals)
Itemp_high=sum(L2_data_integrals);

L1_errors=[L1_errors; abs(Itemp_high-Itemp_low)];








function Ihigh=Ihigh_computation(L2_data)

%--------------------------------------------------------------------------
% Object:
% Approximation of the integral in view of approximations in LIST 2.
%--------------------------------------------------------------------------
% Input:
% L2_data: struct defining LIST 2 sph. triangles vertices, areas and
%          integral approximations.
%--------------------------------------------------------------------------
% Output:
% Ihigh:   approximation of integral of "f" over the initial spherical
%          triangle.
%--------------------------------------------------------------------------

M=length(L2_data);
Ihigh=0;
for j=1:M
    L2_data_temp=L2_data{j};
    L2_data_temp_sons_integral=zeros(4,1);
    for jj=1:4
        L2_data_temp_sons_integral(jj)=L2_data_temp{jj}.integral;
    end
    Ihigh=Ihigh+sum(L2_data_temp_sons_integral);
end











function L2_data=generate_sphtri_sons(vertices,f)

%--------------------------------------------------------------------------
% 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:
% L2_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];

L2_data{1}=make_L2_data(vertices,f);
L2_data{2}=make_L2_data([OAB_mid; OB; OBC_mid],f);
L2_data{3}=make_L2_data([OBC_mid; OC; OAC_mid],f);
L2_data{4}=make_L2_data([OAB_mid; OBC_mid; OAC_mid],f);









function L2_dataL=make_L2_data(vertices,f)

%--------------------------------------------------------------------------
% 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:
% L2_dataL : determine from the spherical triangle defined by vertices, a
%            struct data, with form:
%            * L2_dataL.vertices,
%            * L2_dataL.area,
%            * L2_dataL.integral.
%--------------------------------------------------------------------------

% vertices
L2_dataL.vertices=vertices;
P1=(vertices(1,:))'; P2=(vertices(2,:))'; P3=(vertices(3,:))';
% xyzw = sphtriquad(4,5,P1,P2,P3,0);
xyzw= cub_sphtri(4,P1,P2,P3);
L2_dataL.integral=(xyzw(:,4))'*feval(f,xyzw(:,1),xyzw(:,2),xyzw(:,3));










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


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










function xyzw = cub_sphtri(n,P1,P2,P3)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% The routine computes a cubature formula "xyzw" with (numerical) algebraic
% degree of precision "n" on the spherical triangle with vertices "P1",
% "P2", "P3", by a 2D formula of degree of precision "m+n" on the xy
% projection of the spherical triangle rotated to put the barycenter at the
% north pole.
%
% The parameter "m" depends on the size of the circumradius and it is "2"
% for "small" "r", becoming bigger as "r" approaches "1".
%
% IMPORTANT: The Chebfun toolbox must be installed.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% n: algebraic degree of precision of the rule;
% P1,P2,P3: column arrays of the spherical triangle vertices coords
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% xywz : 4-column array of nodes cartesian coords and weights
%--------------------------------------------------------------------------
% Note:
%--------------------------------------------------------------------------
% The code has a warning at line 126, difficult to skip due to a matrix
% whose size cannot easily be forecasted before the process.
%--------------------------------------------------------------------------
% Routines called:
%--------------------------------------------------------------------------
% 1. compute_m (attached)
% 2. cub_circsect (attached)
%
% The routine "cub_circsect" requires "quad_trig", "r_jacobi" and "gauss"
% that are attached below as well as the "quad_trig" subroutines.
%--------------------------------------------------------------------------
%% 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 2 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, Marco Vianello.
%%
%% Date: January 01, 2021
%% Update: December 29, 2025
%--------------------------------------------------------------------------


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

% working with vertices as column vectors
if size(P1,1) == 1, P1=P1'; end
if size(P2,1) == 1, P2=P2'; end
if size(P3,1) == 1, P3=P3'; end

% ------------------------------- main code -------------------------------
% We scale the problem in the unit-sphere with center [0 0 0].
% In other words, we make up a cubature rule 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).

% ... barycenter ...
RR=norm(P1);
vert0=[P1 P2 P3]/RR;
CC=1/3*(P1+P2+P3); CC=CC/norm(CC);

% .................... rotation matrix to the north pole ..................

[az,el] = 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 ..............

vert1=R*vert0;

% ...... computing "m" needed in determining the degree of precision ......
m=compute_m(vert1');

% .................... determining quadrature rule ........................

xyzw=[]; % nodes on the sphere
vert1=[vert1 vert1(:,1)];

for i=1:3

    % affine transformation matrix
    P=vert1(:,i); Q=vert1(:,i+1);
    om=acos(P'*Q);
    xi=cos(om/2); eta=sin(om/2);
    M=[xi eta 0 0; 0 0 xi eta; xi -eta 0 0; 0 0 xi -eta];
    h=[P(1); P(2); Q(1); Q(2)];
    u=M\h;
    T=[u(1) u(2); u(3) u(4)];

    % nodes and weights for the xy-projection of the rotated spher.triangle
    nw=cub_circsect(m+n,om/2,0,1); nod2=T*nw(:,1:2)';

    % spherical triangle nodes and weights on North Pole spher.triangle
    x=nod2(1,:); y=nod2(2,:); z=sqrt(1-x.^2-y.^2);

    weights=abs(det(T))*nw(:,3)./z';

    % inverse rotation to original spherical triangle
    nodes=invR*[x; y; z];

    % cubature rule update
    xyzw=[xyzw;[nodes' weights]];

end

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






function m=compute_m(vert)

%--------------------------------------------------------------------------
% Object:
% It computes the "m" positive integer value depending from the sph. tri.
% with vertices "vert", so that a WAM over the sph. triangle with degree
% equal to "n" can be obtained via WAM of degree "m+n" on its projection on
% the xy-plane.
%--------------------------------------------------------------------------
% Input:
% vert: points defining the sph. triangle (the k-th point is described
%    by the k-th row.
%--------------------------------------------------------------------------
% Output:
% m: positive integer value depending from the sph.triangle with vertices
% "vert", so that a WAM over the sph. triangle with degree equal to "n" can
% be obtained via cubature of degree "m+n" on its projection on the
% xy-plane.
%--------------------------------------------------------------------------
% Data:
% First version: 13/01/2021 by A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------

r=sqrt( (vert(:,1)).^2 + (vert(:,2)).^2 );
r0=max(r);

intvf=[0,r0]; 
F=@(r) sqrt(1-r); 
f=chebfun(F,intvf);
m=max(1,(length(f)-1));










function xyw = cub_circsect(n,omega,r1,r2)

%--------------------------------------------------------------------------
% Object:
% The routine computes the nodes and weights of a product gaussian
% formula on a circular annular sector centered at the origin
% with angles in [-omega,omega]
%--------------------------------------------------------------------------
% Input:
% n: algebraic degree of exactness
% omega: half-length of the angular interval, 0<omega<=pi
% r1,r2: internal and external radius, 0<=r1<r2
%--------------------------------------------------------------------------
% Output:
% xyw: (ceil((n+2)/2) x (n+1)) x 3 array of (xnodes,ynodes,weights)
%--------------------------------------------------------------------------
% Required routines:
% 1. r_jacobi.m (www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html)
% 2. gauss.m (www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html)
% 3. quad_trig.m
%--------------------------------------------------------------------------
% Written by Gaspare Da Fies and Marco Vianello, University of Padova
% Date: November 8, 2011.
% Last update: January 4. 2020.
%--------------------------------------------------------------------------

% trigonometric gaussian formula on the arc
tw=quad_trig(n,-omega,omega);

% algebraic gaussian formula on the radial segments
ab=r_jacobi(ceil((n+2)/2),0,0);
xw=gauss(ceil((n+2)/2),ab);
xw(:,1)=xw(:,1)*(r2-r1)/2+(r2+r1)/2;
xw(:,2)=xw(:,2)*(r2-r1)/2;

% creating the polar grid
[r,theta]=meshgrid(xw(:,1),tw(:,1));
[w1,w2]=meshgrid(xw(:,2),tw(:,2));

% nodal cartesian coordinates and weights
xyw(:,1)=r(:).*cos(theta(:));
xyw(:,2)=r(:).*sin(theta(:));
xyw(:,3)=r(:).*w1(:).*w2(:);










function tw=quad_trig(n,alpha,beta)

%--------------------------------------------------------------------------
% Object:
% Computation of trigonometric gaussian rules on the unit arc from "alpha"
% to "beta".
%--------------------------------------------------------------------------
% Inputs:
% n    : the rule computes computes the n+1 angles and weights of a
%        trigonometric gaussian quadrature formula on [alpha,beta], with
%                          0 < beta-alpha <= 2*pi;
% alpha, beta: arc angles for trigonometric gaussian rules on the unit arc
%       from "alpha" to "beta".
%--------------------------------------------------------------------------
% Outputs
% tw   : (n+1) x 2 matrix, where the first column contains the nodes, while
%        the second one contains the weights.
%--------------------------------------------------------------------------
% First version: May 18, 2013 by G. Da Fies, A. Sommariva, M. Vianello.
% Successive versions have been made by G. Meurant, A. Sommariva and
% M. Vianello.
% Last release: January 04, 2020.
%--------------------------------------------------------------------------

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

if nargin < 1, n=10; end
if nargin < 2, beta=pi; alpha=-beta; end
if nargin < 3, beta=alpha; alpha=-beta; end

n=n+1;
omega=(beta-alpha)/2;
ab = r_subchebyshev(n,omega);
xw_symm_eigw = SymmMw(n,ab);
tw=quad_trig_conversion(xw_symm_eigw,omega);
tw(:,1)=tw(:,1)+(beta+alpha)/2;









function ab=r_subchebyshev(n,omega)

%--------------------------------------------------------------------------
% Object:
% Recurrence coeffs for the monic OPS w.r.t. the weight function
%         w(x)=2*sin(omega/2)/sqrt(1-sin^2(omega/2)*x^2)
% by the modified Chebyshev algorithm.
% The reference angle of the rule is [-omega,omega].
%--------------------------------------------------------------------------
% Inputs:
% n     : number of points.
% omega : arc angle.
%--------------------------------------------------------------------------
% Output:
% ab   : three terms recursion
%--------------------------------------------------------------------------
% First version: May 18, 2013 by G. Da Fies, A. Sommariva, M. Vianello.
% Successive versions have been made by G. Meurant, A. Sommariva and
% M. Vianello.
% Last release: January 04, 2020.
%--------------------------------------------------------------------------

N=n; n=n-1;

% modified Chebyshev moments by recurrence
if rem(N,2) == 1
    NN=N+1; nn=n+1;
else
    NN=N; nn=n;
end
mom=fast_moments_computation(omega,2*nn+1);

% recurrence coeffs of the monic Chebyshev polynomials
abm(:,1)=zeros(2*nn+1,1);
abm(:,2)=0.25*ones(2*nn+1,1); abm(1,2)=pi; abm(2,2)=0.5;

% recurrence coeffs for the monic OPS w.r.t. the weight function
ab = fast_chebyshev(NN,mom,abm);









function x = tridisolve(a,b,c,d)

%--------------------------------------------------------------------------
% Object:
% Solution of tridiagonal system of equations.
%--------------------------------------------------------------------------
% From Cleve Moler's Matlab suite
% http://www.mathworks.it/moler/ncmfilelist.html
%--------------------------------------------------------------------------
% x = TRIDISOLVE(a,b,c,d) solves the system of linear equations
%
%     b(1)*x(1) + c(1)*x(2) = d(1),
%     a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,
%     a(n-1)*x(n-1) + b(n)*x(n) = d(n).
%
% The algorithm does not use pivoting, so the results might be inaccurate
% if abs(b) is much smaller than abs(a)+abs(c).
%
% More robust, but slower, alternatives with pivoting are:
%     x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)
%     x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)
%--------------------------------------------------------------------------
% Optimized version by G. Meurant.
%--------------------------------------------------------------------------

x = d;
n = length(x);
bi = zeros(n,1);

for j = 1:n-1
    bi(j) = 1 / b(j);
    mu = a(j) * bi(j);
    b(j+1) = b(j+1) - mu * c(j);
    x(j+1) = x(j+1) - mu * x(j);
end

x(n) = x(n) / b(n);
for j = n-1:-1:1
    x(j) = (x(j) - c(j) * x(j+1)) * bi(j);
end






function ab=fast_chebyshev(N,mom,abm)

%--------------------------------------------------------------------------
% Object:
% Modified Chebyshev algorithm, that works only for the subperiodic weight
% function.
%--------------------------------------------------------------------------
% From Gautschi's code (simplified)
% Mar 2012
%--------------------------------------------------------------------------
% Optimized version by G. Meurant.
%--------------------------------------------------------------------------

ab = zeros(N,2);
sig = zeros(N+1,2*N);

ab(1,2) = mom(1);

sig(1,1:2*N) = 0;
sig(2,:) = mom(1:2*N);

for n = 3:N+1
    for m = n-1:2*N-n+2
        sig(n,m) = sig(n-1,m+1) + abm(m,2) * sig(n-1,m-1) - ...
            ab(n-2,2) * sig(n-2,m);
    end

    ab(n-1,2) = sig(n,n-1) / sig(n-1,n-2);
end









function mom=fast_moments_computation(omega,n)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Inputs:
%--------------------------------------------------------------------------
% Outputs:
%--------------------------------------------------------------------------
% Authors G. Meurant and A. Sommariva
% June 2012
%--------------------------------------------------------------------------

mom=zeros(1,n+1);
mom(1)=2*omega; % FIRST MOMENT.

if(n>=2)

    if(omega<=1/4*pi)
        l=10;
    elseif(omega<=1/2*pi)
        l=20;
    elseif(omega<=3/4*pi)
        l=40;
    else
        if omega == pi
            l=2*ceil(10*pi);
        else
            l=2*ceil(10*pi/(pi-omega));
        end
    end


    temp=(2:2:n+2*l-2); % AUXILIAR VECTORS.
    temp2=temp.^2-1;

    dl=1/4 -1./(4*(temp-1)); % DIAGONALS.
    dc=1/2 -1/sin(omega/2)^2 -1./(2*temp2);
    du=1/4 +1./(4*(temp+1));

    d=4*cos(omega/2)/sin(omega/2)./temp2'; % COMPUTING KNOWN TERM.
    d(end)=d(end);                         % PUT LAST MOMENT NULL

    z=tridisolve(dl(2:end),dc,du(1:end-1),d); % SOLVE SYSTEM.
    mom(3:2:n+1)=z(1:floor(n/2)); % SET ODD MOMENTS.

end

mom=mom';

normalized = 0;

if normalized == 0
    M=length(mom);
    kk=2.^(-((1:2:M)-2))'; kk(1)=1;
    v=ones(M,1);
    v(1:2:M)=kk;
    mom=v.*mom;
end









function xw=SymmMw(N,ab)

%--------------------------------------------------------------------------
% Object:
% Computation of the nodes and weights for a symmetric weight function
% this version uses the reduced matrix and eig and computation of weights
% with the 3-term recurrence.
%--------------------------------------------------------------------------
% Input:
% N : cardinality of the rule
% ab: 3-term recurrence for the orthogonal polynomials same as in OPQ
%     ab(1,2) is the 0th moment.
%--------------------------------------------------------------------------
% Output:
% xw : xw(:,1) nodes, xw(:,2) weights of the quadrature rule
%--------------------------------------------------------------------------
% Reference paper:
% Fast variants of the Golub and Welsch algorithm for symmetric
% weight functions by G. Meurant and A. Sommariva (2012)
%--------------------------------------------------------------------------
% Data:
% Written by G. Meurant and A. Sommariva on June 2012
%--------------------------------------------------------------------------

N0 = size(ab,1);
if N0 < N
    error('SymmMw: input array ab is too short')
end

na = norm(ab(:,1));
if na > 0
    error('SymmMw: the weight function must be symmetric')
end

% computation of the reduced matrix in vectors (a,b)

if mod(N,2) == 0
    even = 1;
    Nc = N / 2;
else
    even = 0;
    Nc = fix(N / 2) +1;
end


absd = ab(:,2);
absq = sqrt(absd);

a = zeros(1,Nc);
b = a;

switch even
    case 1
        % N even
        a(1) = absd(2);
        b(1) = absq(2) * absq(3);

        k = 2:Nc-1;
        a(k) = absd(2*k-1) + absd(2*k);
        b(k) = absq(2*k) .* absq(2*k+1);
        a(Nc) = absd(N) + absd(N-1);
        start = 1;

        J = diag(a) + diag(b(1:Nc-1),1) + diag(b(1:Nc-1),-1);
        t = sort(eig(J));
        w = weights_3t(t',a,b);
        % w are the squares of the first components
        w = w' / 2;
    case 0
        % N odd
        a(1) = absd(2);
        b(1) = absq(2) * absq(3);

        k = 2:Nc-1;
        a(k) = absd(2*k-1) + absd(2*k);
        b(k) = absq(2*k) .* absq(2*k+1);
        a(Nc) = absd(N);
        start = 2;

        % the first node must be zero
        J = diag(a) + diag(b(1:Nc-1),1) + diag(b(1:Nc-1),-1);
        t = sort(eig(J));
        t(1) = 0;
        w = weights_3t(t',a,b);
        w = [w(1); w(2:end)' / 2];
    otherwise
        error('this is not possible')
end

xwp = sqrt(t);

xw(:,1) = [-xwp(end:-1:start,1); xwp];
xw(:,2) = ab(1,2) * ([w(end:-1:start); w]);









function tw=quad_trig_conversion(xw,omega)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Inputs:
%--------------------------------------------------------------------------
% Outputs:
%--------------------------------------------------------------------------
% Authors G. Meurant and A. Sommariva
% June 2012
%--------------------------------------------------------------------------

tw(:,1)=2*asin(sin(omega/2)*xw(:,1));
tw(:,2)=xw(:,2);









function w=weights_3t(t,a,b)

%--------------------------------------------------------------------------
% Object:
% Squares of the 1st components of eigenvectors from the 3-term
% recurrence relation of the orthogonal polynomials
%--------------------------------------------------------------------------
% Inputs:
% t: nodes
% a,b: coefficients of the 3-term recurrence
%--------------------------------------------------------------------------
% Outputs
% w: squares of the first components of the eigenvectors
%--------------------------------------------------------------------------
% Authors G. Meurant and A. Sommariva
% June 2012
%--------------------------------------------------------------------------

N = length(t);

P = zeros(N,N);
P(1,:) = ones(1,N);
P(2,:) = (t - a(1)) / b(1);

for k = 3:N
    k1 = k - 1;
    k2 = k - 2;
    P(k,:) = ((t - a(k1)) .* P(k1,:) - b(k2) * P(k2,:)) / b(k1);
end

P2 = P .* P;
w = 1 ./ sum(P2);










% GAUSS Gauss quadrature rule.
%
%    Given a weight function w encoded by the nx2 array ab of the
%    first n recurrence coefficients for the associated orthogonal
%    polynomials, the first column of ab containing the n alpha-
%    coefficients and the second column the n beta-coefficients,
%    the call xw=GAUSS(n,ab) generates the nodes and weights xw of
%    the n-point Gauss quadrature rule for the weight function w.
%    The nodes, in increasing order, are stored in the first
%    column, the n corresponding weights in the second column, of
%    the nx2 array xw.
%
function xw=gauss(N,ab)
N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
    J(n,n-1)=sqrt(ab(n,2));
    J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];









function ab=r_jacobi(N,a,b)

nu=(b-a)/(a+b+2);
mu=2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2);
if N==1
    ab=[nu mu]; return
end

N=N-1;
n=1:N;
nab=2*n+a+b;
nuadd=(b^2-a^2)*ones(1,N)./(nab.*(nab+2));
A=[nu nuadd];
n=2:N;
nab=nab(n);
B1=4*(a+1)*(b+1)/((a+b+2)^2*(a+b+3));
B=4*(n+a).*(n+b).*n.*(n+a+b)./((nab.^2).*(nab+1).*(nab-1));
abadd=[mu; B1; B'];
ab=[A' abadd];




















