
function [thetaw,xw]=cub_sphrect(n,sphrect_intv,method)


%--------------------------------------------------------------------------
% Object:
%
% Cubature on spherical quadrangles determined by "sphrect_intv".
% The cubature rule has algebraic degree of precision "n".
%--------------------------------------------------------------------------
% Input:
%
% n: degree of algebraic polynomial exactness of the rule.
%
% 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]).
%
% method: 0: Gaussian rule, 1 Fejer I, 2 Fejer II, 3 Clenshaw-Curtis.
%
%--------------------------------------------------------------------------
% 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.
%
% xw:      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. "quad_trig.m"(EXTERNAL TO THIS FILE),
% 2. "quad_trig_abssin.m" (ATTACHED TO THIS FILE),
% 3. "quad_trig_cc.m" (EXTERNAL TO THIS FILE),
% 4. "circle_trigquad" (ATTACHED TO THIS FILE).
%--------------------------------------------------------------------------
%% 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@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: May 17, 2013
%% Last Update: December 21, 2020.
%--------------------------------------------------------------------------

 % ......................... troubleshooting ..............................
 
 if nargin < 3, method=0; end
 
 
 
 
% ......................... 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=[]; xw=[];
else
    switch method
        
        % ....................... quad_trig .......................
        
        case 0
            
            % .................... 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,sphrect_intv(1,1),...
                    sphrect_intv(1,2));
                yy(2,1:n+1,:)=quad_trig_abssin(n,sphrect_intv(2,2));
                
                % In the case of a cap we will consider the volume
                % element directly in quad_trig_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,:)=quad_trig(n,sphrect_intv(1,1),...
                    sphrect_intv(1,2));
                yy(2,1:n+2,:)=quad_trig(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
            
        otherwise
            
            % ................ Weighted Fejer / Clenshaw Curtis ...........
            
            [x1,w1]=quad_trig_cc(2*n+1,sphrect_intv(1,1),...
                sphrect_intv(1,2),method);
            [x2,w2]=quad_trig_cc(2*n+2,sphrect_intv(2,1),...
                sphrect_intv(2,2),method);
            [A,B]=meshgrid(x1,x2);
            [C,D]=meshgrid(w1,w2);
            A=A(:); B=B(:); C=C(:); D=D(:);
            ww=C.*D.*sin(B);
            thetaw=[A B ww];
            
    end
    
end

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

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

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

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









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

function tw=circle_trigquad(n,alpha,beta)

if nargin == 1
    alpha=0;
    beta=2*pi;
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 tw=quad_trig_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=quad_trig_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=quad_trig_fast_sin(n,omega);
tw(:,1)=0.5*(alpha+beta)+tw(:,1);




%--------------------------------------------------------------------------
% ATTACHED ROUTINES.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% quad_trig_fast_sin
%--------------------------------------------------------------------------

function tw=quad_trig_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);
xw=SymmMw(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






%--------------------------------------------------------------------------
% SymmMw
%--------------------------------------------------------------------------

function xw=SymmMw(N,ab)
%SYMMMw 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

%
% see: Fast variants of the Golub and Welsch algorithm for symmetric
% weight functions by G. Meurant and A. Sommariva (2012)

% 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
% nloop: number of iterations in QR

%
% Authors G. Meurant and A. Sommariva
% 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]);






%--------------------------------------------------------------------------
% weights_3t
%--------------------------------------------------------------------------

function w=weights_3t(t,a,b);
%WEIGHTS_3T squares of the 1st components of eigenvectors from the 3-term
% recurrence relation of the orthogonal polynomials
%

% Input
% t: nodes
% a,b coefficients of the 3-term recurrence
%
% Ouput
% 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);






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








