

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. trigauss.m 
%
% All these procedures are attached below. 
%--------------------------------------------------------------------------
% Reference paper
%--------------------------------------------------------------------------
% Algebraic cubature by linear blending of elliptical arcs,
% Applied Numerical Mathematics, 74 (2013), 49-61. (Open Access)
% G. Da Fies, A. Sommariva and M. Vianello.
%--------------------------------------------------------------------------
% Written by Gaspare Da Fies, Alvise Sommariva and Marco Vianello.
%
% Date: November 8, 2011.
% Last update: January 4. 2026.
%--------------------------------------------------------------------------
% Copyright
%--------------------------------------------------------------------------
%% Copyright (C) 2013-
%% Gaspare Da Fies, 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:
%% Gaspare Da Fies, Alvise Sommariva, Marco Vianello.
%%
%% Date: MAY 18, 2013
%% Modified on: JANUARY 4, 2026
%--------------------------------------------------------------------------

% trigonometric gaussian formula on the arc
tw=trigauss(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(:);









%--------------------------------------------------------------------------
% Subroutines
%--------------------------------------------------------------------------

function [tw,xw,ab]=trigauss(tde,alpha,beta)

%--------------------------------------------------------------------------
% Object
%--------------------------------------------------------------------------
% Computation of trigonometric gaussian rules on the unit arc 
% "[alpha,beta]" with trigonometric degree of exactness "tde".
%--------------------------------------------------------------------------
%            Original note by G. Da Fies and M. Vianello:
%--------------------------------------------------------------------------
% From the original code by G. Da Fies and M. Vianello:
%
% computes the n+1 angles and weights of a trigonometric gaussian
% quadrature formula on [alpha,beta], 0<beta-alpha<=pi
%
% uses the routines chebyshev.m, gauss.m from
% www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html
% we suggest to put the following statements
% ab = zeros(N,2); sig = zeros(N+1,2*N);
% at the beginning of the body of chebyshev.m to speed-up execution
%
% the formula integrates the canonical trigonometric basis with accuracy
% from about 10^(-15) (small omega) to about 10^(-13) (omega-->pi)
% up to n=300
%--------------------------------------------------------------------------
% Inputs
%--------------------------------------------------------------------------
% n    : the rule computes computes the "tde+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". Note that "0<beta-alpha<=2*pi"
%--------------------------------------------------------------------------
% Outputs
%--------------------------------------------------------------------------
% tw   : (tde+1) x 2 matrix, where  
%        a) "tw(:,1)" determines the nodes of the trig. rule, 
%        b) "tw(:,2)" determines the weights of the trig. rule
% xw   : (tde+1) x 2 matrix, where  
%        a) "xw(:,1)" determines the nodes of the associated alg. rule, 
%        b) "xw(:,2)" determines the weights of the associated alg. rule,
% ab   : (tde+1) x 2 matrix, storing coefficients of the Jacobi matrix of
%        the associated alg. rule.
%--------------------------------------------------------------------------
% Versions of this code
%--------------------------------------------------------------------------
% 1. First version: November 8, 2011
%    By G. Da Fies, and M. Vianello.
% 
% 2. Modified on: December 31, 2025 (by A. Sommariva).
%--------------------------------------------------------------------------
% Example
%--------------------------------------------------------------------------
% >> tde=5; alpha=pi/4; beta=pi/2;
% >> [tw,xw,ab]=trigauss_classic_2011(tde,alpha,beta);
% >> tw
% 
% tw =
% 
%     0.8122    0.0678
%     0.9192    0.1419
%     1.0848    0.1830
%     1.2714    0.1830
%     1.4370    0.1419
%     1.5440    0.0678
% 
% >> xw
% 
% xw =
% 
%    -0.9327    0.0678
%    -0.6618    0.1419
%    -0.2390    0.1830
%     0.2390    0.1830
%     0.6618    0.1419
%     0.9327    0.0678
% 
% >> ab
% 
% ab =
% 
%          0    0.7854
%          0    0.3351
%          0    0.2663
%          0    0.2571
%          0    0.2539
%          0    0.2525
% 
% >> 
%--------------------------------------------------------------------------
%% Copyright (C) 2013-
%% Gaspare Da Fies, 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:
%% Gaspare Da Fies, Alvise Sommariva, Marco Vianello.
%%
%% Date: MAY 18, 2013
%% Date: DECEMBER 31, 2025
%--------------------------------------------------------------------------

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

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

% n: cardinality of the rule
n=tde+1;

% half-length of the angular interval
omega=(beta-alpha)/2;

if omega == pi

    ab=[];
    tw=circle_trigquad(n,alpha);
    xw=[];

else
    
    ab=r_trigauss(n-1,omega);
    
    % Gaussian formula for the weight function above
    xw=gauss(n,ab);

    % angles and weights for the trigonometric gaussian formula
    tw(:,1)=2*asin(sin(omega/2)*xw(:,1))+(beta+alpha)/2;
    tw(:,2)=xw(:,2);
    
end






%--------------------------------------------------------------------------
% Subroutines
%--------------------------------------------------------------------------

function ab=r_trigauss(n,omega)

% modified Chebyshev moments by recurrence
z(1)=2*omega;

z(n+1)=quadgk(@(t)cos(2*n*acos(sin(t/2)/sin(omega/2))),...
    -omega,omega,'MaxIntervalCount',50000);

temp=(2:2:2*n-1);
dl=1/4-1./(4*(temp-1));

dc=1/2-1/sin(omega/2)^2-1./(2*(temp.^2-1));
du=1/4+1./(4*(temp+1));
d=4*cos(omega/2)/sin(omega/2)./(temp.^2-1)';
d(n-1)=d(n-1)-du(n-1)*z(n+1);
z(2:n)=tridisolve(dl(2:n-1),dc(1:n-1),du(1:n-2),d(1:n-1));
mom=zeros(1,2*n+2);
mom(1:2:2*n+1)=z(1:n+1);

% normalization of the moments (monic polynomials)
k=(3:length(mom));
mom(3:end)=exp((2-k)*log(2)).*mom(3:end);

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

% 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
ab=chebyshev(n+1,mom,abm);









function x = tridisolve(a,b,c,d)
%   TRIDISOLVE  Solve 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)

x = d;
n = length(x);
for j = 1:n-1
    mu = a(j)/b(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))/b(j);
end










function [ab,normsq]=chebyshev(N,mom,abm)

% CHEBYSHEV Modified Chebyshev algorithm.
%
%    Given a weight function w encoded by its first 2n modified
%    moments, stored in the (row) vector mom, relative to monic 
%    polynomials defined by the (2n-1)x2 array abm of their
%    recurrence coefficients, [ab,normsq]=CHEBYSHEV(n,mom,abm)
%    generates the array ab of the first n recurrence coefficients
%    of the orthogonal polynomials for the weight function w, and 
%    the vector normsq of their squared norms. The n alpha-
%    coefficients are stored in the first column, the n beta-
%    coefficients in the second column, of the nx2 array ab. The
%    call [ab,normsq]=CHEBYSHEV(n,mom) does the same, but using the 
%    classical Chebyshev algorithm. If n is larger than the sizes
%    of mom and abm warrant, then n is reduced accordingly.
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.

if N<=0, error('N out of range'), end
if N>size(mom,2)/2, N=size(mom,2)/2; end
if nargin<3, abm=zeros(2*N-1,2); end
if N>(size(abm,1)+1)/2, N=(size(abm,1)+1)/2; end
ab(1,1)=abm(1,1)+mom(2)/mom(1); ab(1,2)=mom(1);
if N==1, normsq(1)=mom(1); return, end
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)-(ab(n-2,1)-abm(m,1))*sig(n-1,m) ...
      -ab(n-2,2)*sig(n-2,m)+abm(m,2)*sig(n-1,m-1);
  end
  ab(n-1,1)=abm(n-1,1)+sig(n,n)/sig(n,n-1)-sig(n-1,n-1)/ ...
    sig(n-1,n-2);
  ab(n-1,2)=sig(n,n-1)/sig(n-1,n-2);
end

normsq=zeros(N,1);

for n=1:N
    normsq(n)=sig(n+1,n); 
end
normsq=normsq';









function xw=gauss(N,ab)

% 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.
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.

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 tw=circle_trigquad(n,alpha)

if nargin < 2
    alpha=0;
end

N=n+1;
w=(2*pi/N)*ones(N,1);
t=linspace(pi/N,2*pi-pi/N,N); t=alpha+t';

tw=[t w];









function ab=r_jacobi(N,a,b)

% R_JACOBI Recurrence coefficients for monic Jacobi polynomials.
%
%    ab=R_JACOBI(n,a,b) generates the first n recurrence
%    coefficients for monic Jacobi polynomials with parameters
%    a and b. These are orthogonal on [-1,1] relative to the
%    weight function w(t)=(1-t)^a(1+t)^b. The n alpha-coefficients
%    are stored in the first column, the n beta-coefficients in
%    the second column, of the nx2 array ab. The call ab=
%    R_JACOBI(n,a) is the same as ab=R_JACOBI(n,a,a) and
%    ab=R_JACOBI(n) the same as ab=R_JACOBI(n,0,0).
%
%    Supplied by Dirk Laurie, 6-22-1998; edited by Walter
%    Gautschi, 4-4-2002.

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];