function [thetaw,xyzw]=cub_sphrect(sphrect_intv,n)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Cubature on spherical quadrangles of longitude-latitude type determined
% by "sphrect_intv".
%
% The cubature rule has algebraic degree of precision "ade".
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% sphrect_intv: matrix 2 x 2;
%     Horizontally, "sphrect_intv" must denote the start point and the end
%     point for each range.
%     Notice that
%     * the first interval must be contained in [0,pi],
%     * the second interval must be contained in [0,2*pi]).
%
%     Notice that each point will be represented as
%
%           x_1=sin(theta_1)*cos(theta_2),
%           x_2=sin(theta_1)*sin(theta_2),
%           x_3=cos(theta_1)
%
%     with "theta1" in "[0,pi]" and "theta2" in "[0,2*pi]"
%
% n: degree of algebraic polynomial exactness of the rule.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% thetaw: For each row it will denote the spherical coordinates of
%          nodes, i.e.
%
%           x_1=sin(theta_1)*cos(theta_2),
%           x_2=sin(theta_1)*sin(theta_2),
%           x_3=cos(theta_1).
%
%          (with theta_1 in [0,pi]), theta_2 in [0,2*pi] and in the last
%          column it will store the weights for each node.
%
% xyzw:      For each row it will denote the cartesian coordinates of
%          nodes, in the last column it will store the weights for each
%          node.
%--------------------------------------------------------------------------
% Required routines:
%
% This program will use:
%
% 1. circle_trigquad (attached below)
% 2. trigauss_abssin
% 3. trigauss
%--------------------------------------------------------------------------
%% Copyright (C) 2013. Mariano Gentile, 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:
%%          Mariano Gentile  <gentimar@hotmail.com>
%%          Alvise Sommariva <alvise@math.unipd.it>
%%          Marco Vianello   <marcov@math.unipd.it>
%%
%% Date       : May 17, 2013
%% Last Update: January 3, 2026
%--------------------------------------------------------------------------


% ......................... Main code below .-.............................

% Using alternative notation.
sphrect_intv=[sphrect_intv(2,:); sphrect_intv(1,:)];

d=size(sphrect_intv,1);
if d~=2
    fprintf('WARNING: Use cubature_d_sphere\n'); thetaw=[]; xyzw=[];
else
    % .................... SPHERICAL CAP ..........................

    if (sphrect_intv(1,2)-sphrect_intv(1,1)==2*pi) && ...
            (sphrect_intv(2,1)==0)
        if mod(n,2)==1, nn=n+1; else, nn=n+2; end
        yy=ones(d,nn,2);
        % We need even number of nodes in [0,2*pi] to have a
        % central simmetry.
        yy(1,1:nn,:)=circle_trigquad(nn-1);
        yy(2,1:n+1,:)=trigauss_abssin(n,sphrect_intv(2,2));

        % In the case of a cap we will consider the volume
        % element directly in trigauss_cap.
        [A,B]=meshgrid(yy(1,1:nn,1),yy(2,1:floor((n+1)/2),1));
        yys(:,1)=A(:);
        yys(:,2)=B(:);
        [C,D]=meshgrid(yy(1,1:nn,2),yy(2,1:floor((n+1)/2),2));
        yys(:,3)=C(:).*D(:);
        % Manually add the North Pole.
        if mod(n,2)==0
            yys(end+1,2)=0; yys(end,3)=yy(2,floor((n+1)/2)+1,2)*pi;
        end

        thetaw=yys;

        % ................. SPHERICAL QUADRANGLE ..................

    else
        yy=ones(d,n+2,2);
        yy(1,1:n+1,:)=trigauss(n,sphrect_intv(1,1),...
            sphrect_intv(1,2));
        yy(2,1:n+2,:)=trigauss(n+1,sphrect_intv(2,1),...
            sphrect_intv(2,2));
        [A,B]=meshgrid(yy(1,1:n+1,1),yy(2,1:n+2,1));
        yys(:,1)=A(:);
        yys(:,2)=B(:);
        [C,D]=meshgrid(yy(1,1:n+1,2),yy(2,1:n+2,2));
        yys(:,3)=C(:).*D(:).*sin(yys(:,2));
        thetaw=yys;

    end

end


% ................... Pointset in cartesian coordinates ...................

if nargout > 1
    xyzw=ones(size(thetaw,1),3);
    xyzw(:,1)=cos(thetaw(:,1)).*sin(thetaw(:,2)); % NODES.
    xyzw(:,2)=sin(thetaw(:,1)).*sin(thetaw(:,2));
    xyzw(:,3)=cos(thetaw(:,2));
    xyzw(:,4)=thetaw(:,3); % WEIGHTS.
end

% ................... Pointset in spherical coordinates ...................

thetaw=[thetaw(:,2) thetaw(:,1) thetaw(:,3)];









%--------------------------------------------------------------------------
% ADDITIONAL ROUTINES.
%--------------------------------------------------------------------------

function tw=circle_trigquad(n)

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

tw=[t w];









function tw=trigauss_abssin(n,omega)

%--------------------------------------------------------------------------
% PURPOSE.
%--------------------------------------------------------------------------
% This routine computes the n+1 angles and weights of a trigonometric
% gaussian quadrature formula on [-omega,omega], 0<omega<=pi,
% with respect to the weight function w(theta)=abs(sin(theta)).
%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% n: trigonometric degree of exactness
% omega: determines the symmetric angular interval [-omega,omega]
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% tw: (n+1) x 2 array of (angles,weights)
%--------------------------------------------------------------------------
% EXAMPLE.
%--------------------------------------------------------------------------
% >> tw=trigauss_abssin(5,pi/4)
%
% tw =
%
%    -0.7376    0.0814
%    -0.5480    0.1302
%    -0.2577    0.0814
%     0.2577    0.0814
%     0.5480    0.1302
%     0.7376    0.0814
%
% >>
%--------------------------------------------------------------------------
%% Copyright (C) 2013
%% Mariano Gentile, 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:
%% Mariano Gentile, Alvise Sommariva, Marco Vianello.
%%
%% Date: MAY 16, 2013
%--------------------------------------------------------------------------

if nargin < 1, n=10; end
if nargin < 2, omega=pi; end
if omega < 0, omega=-omega; end

beta=omega;
alpha=-omega;

omega=(beta-alpha)/2;
tw=trigauss_fast_sin(n,omega);
tw(:,1)=0.5*(alpha+beta)+tw(:,1);






%--------------------------------------------------------------------------
% trigauss_fast_sin
%--------------------------------------------------------------------------

function tw=trigauss_fast_sin(n,omega)

%--------------------------------------------------------------------------
% PURPOSE.
%-----------------
% This routine computes the n+1 angles and weights of a trigonometric
% gaussian quadrature formula on [-omega,omega], 0<omega<=pi with
% respect to the weight function w(theta)=abs(sin(theta)).
%--------------------------------------------------------------------------
% INPUT.
%-----------------
% n: trigonometric degree of exactness.
% omega: determines the angular interval (-omega,omega).
%--------------------------------------------------------------------------
% OUTPUT.
%-----------------
% tw: (n+1) x 2 array of (angles,weights)
%--------------------------------------------------------------------------
% ROUTINES.
%-----------------
% 1. fast_moment_evaluation (attached to this file),
% 2. fast_chebyshev (attached to this file),
% 3. SymmMw (attached to this file).
%--------------------------------------------------------------------------

% modified Chebyshev moments
mom=fast_moment_evaluation(n,omega);


% 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)=4*(sin(omega/2))^2*abs(x) (it will
% by the modified Chebyshev algorithm
% [ab,normsq]=chebyshev(n+1,mom,abm);
[ab]=fast_chebyshev(n+1,mom,abm);

% Gaussian formula for the weight function above
xw=gauss(n+1,ab);

% angles and weights for the trigonometric gaussian formula w.r.t. the
% weight function w(theta)=abs(sin(theta))
tw(:,1)=2*asin(sin(omega/2)*xw(:,1));
tw(:,2)=xw(:,2);






%--------------------------------------------------------------------------
% fast_moment_evaluation
%--------------------------------------------------------------------------

function mom=fast_moment_evaluation(n,omega)
%--------------------------------------------------------------------------
% PURPOSE.
%--------------------------------------------------------------------------
% This routine computes the moments of Chebyshev polynomials with
% respect to the weight function w(theta)=abs(sin(theta)).
%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% n: trigonometric degree of exactness.
% omega: determines the angular interval (-omega,omega).
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% mom: moments of Chebyshev polynomials w.r.t. w(theta)=abs(sin(theta)).
%--------------------------------------------------------------------------
% ROUTINES.
%--------------------------------------------------------------------------
% 1. expl_mom (attached to this file).
%--------------------------------------------------------------------------
momexpl=expl_mom(n,omega);
mom=zeros(1,2*n+2);
mom(1:2:end)=momexpl;






%--------------------------------------------------------------------------
% fast_chebyshev
%--------------------------------------------------------------------------

function ab=fast_chebyshev(N,mom,abm)
%SUBP_MOD_CHEBYSHEV Modified Chebyshev algorithm
% this works only for the subperiodic weight function
%
% From Gautschi's code (simplified)
% Mar 2012
%

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






%--------------------------------------------------------------------------
% expl_mom
%--------------------------------------------------------------------------

function mom=expl_mom(n,omega)

% Y=cos(2*(k'-1)*acos(x)).*(ones(m,1)*abs(x).*4.*(sin(omega/2))^2);
c=4.*(sin(omega/2))^2;
k=(1:n)';
s=(2*k+1);
ak=(s.*(1-(-1).^(k+1)))./(s.^2-1);
bk=1./(s-1);
it_odd=[1/2; (ak-bk)];
mom=2*[c/2; (c/2)*(it_odd(1:end-1)+it_odd(2:end))];






function tw=trigauss(n,alpha,beta)

%--------------------------------------------------------------------------
% PURPOSE.
%--------------------------------------------------------------------------
% FAST VARIANT OF A CODE BY DA FIES AND VIANELLO FOR THE COMPUTATION OF
% SUB-TRIGONOMETRIC GAUSSIAN RULES.
%--------------------------------------------------------------------------
% INPUTS:
%--------------------------------------------------------------------------
% n    : NUMBER OF POINTS.
% alpha, beta: ARC ANGLES, FOR SUBTRIG. RULE IN [alpha,beta].
%--------------------------------------------------------------------------
% OUTPUTS:
%--------------------------------------------------------------------------
% tw   : QUADRATURE NODES AND WEIGHTS (RESPECTIVELY 1st AND 2nd COLUMN).
%--------------------------------------------------------------------------
% EXAMPLE:
%--------------------------------------------------------------------------
% >> tw=trigauss_2012(5,0,pi/4)
%
% tw =
%
%     0.0268    0.0678
%     0.1338    0.1419
%     0.2994    0.1830
%     0.4860    0.1830
%     0.6516    0.1419
%     0.7586    0.0678
%
% >>
%
%--------------------------------------------------------------------------
%%
%% 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 <alvise@math.unipd.it>
%%          Marco Vianello   <marcov@math.unipd.it>
%%
%% Date     : May 18, 2013
%% Modified : January 3, 2026
%--------------------------------------------------------------------------

if nargin < 1, n=10; end

if (nargin == 2)
    beta=alpha; alpha=-beta;
end

if (nargin == 1)
    beta=pi;
    alpha=-beta;
end

n=n+1;
omega=(beta-alpha)/2;
ab = r_subchebyshev(n,omega);
xw = gauss(n,ab);
tw=trigauss_conversion(xw,omega);
tw(:,1)=tw(:,1)+(beta+alpha)/2;






%--------------------------------------------------------------------------
% r_subchebyshev
%--------------------------------------------------------------------------

function ab=r_subchebyshev(n,omega)

%--------------------------------------------------------------------------
% FAST TRIGAUSS.
% FAST VARIANT OF A CODE BY DA FIES AND VIANELLO.
%--------------------------------------------------------------------------
% INPUTS:
% n    : NUMBER OF POINTS.
% omega: ARC ANGLE.
%
% OUTPUTS:
% ab   : THREE TERMS RECURSION.
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%% Copyright (C) 2007-2012
%% Gaspare Da Fies, Gerard Meurant, 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, Gerard Meurant, Alvise Sommariva, Marco Vianello.
%%
%% Date: June 3, 2012
%--------------------------------------------------------------------------

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
% w(x)=2*sin(omega/2)/sqrt(1-sin^2(omega/2)*x^2)
% by the modified Chebyshev algorithm

% ab = chebyshev(NN+1,mom,abm);
ab = fast_chebyshev(NN,mom,abm);






%--------------------------------------------------------------------------
% tridisolve
%--------------------------------------------------------------------------

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)

% optimized version
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






%--------------------------------------------------------------------------
% fast_moments_computation
%--------------------------------------------------------------------------

function mom=fast_moments_computation(omega,n)

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






%--------------------------------------------------------------------------
% trigauss_conversion
%--------------------------------------------------------------------------

function tw=trigauss_conversion(xw,omega)

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


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