function [thetaw,xw]=cubature_2_subsphere_fast(x,n,fig,method)


%--------------------------------------------------------------------------
% PURPOSE.
%-----------------------------
%
% Cubature on spherical quadrangles determined by "x".
% The cubature rule has algebraic degree of precision n.
%
%--------------------------------------------------------------------------
% INPUT.
%-----------------------------
%
% x = matrix 2 x 2
%     Horizontally, x must denote the start point and the end point for
%     each range.
%     Remember that the first interval must be contained in [0,pi], and
%     the second one in [0,2*pi]).
%
% n = degree of algebraic polynomial exactness of the rule.
%
% method: 0 trigauss, 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. "trigauss_2012.m"(EXTERNAL TO THIS FILE),
% 2. "trigauss_abssin.m" (EXTERNAL TO THIS FILE),
% 3. "trig_cc.m" (EXTERNAL TO THIS FILE),
% 4. "circle_trigquad" (ATTACHED TO THIS FILE).
%
%--------------------------------------------------------------------------
% NOTES.
%-----------------------------
%
%   ORIGINAL CODE BY: Mariano Gentile,
%   REVISED BY: Alvise Sommariva.
%
%   1. This program is optimized for caps in S^2. For general caps in S^d
%      use "cubature_d_sphere.m".
%
%
%%
%% 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
%%
%-------------------------------------------------------------------------


x=[x(2,:); x(1,:)]; % USING GENTILE NOTATION.

if nargin < 3
    fig = 0;
end
if nargin < 4
    method=0;
end

d=size(x,1);
if d~=2
    fprintf('WARNING: Use cubature_d_sphere\n');
    thetaw=[]; xw=[];
else

    switch method
        case 0
            if (x(1,2)-x(1,1)==2*pi) && (x(2,1)==0) % SPHERICAL CAP.
                % fprintf('\n \t CAP');
                if mod(n,2)==1
                    yy=ones(d,n+1,2);
                    % yy(1,1:n+1,:)=trigauss_2012(n,x(1,1),x(1,2));
                    yy(1,1:n+1,:)=circle_trigquad(n,x(1,1),x(1,2));
                    yy(2,1:n+1,:)=trigauss_abssin(n,x(2,2));

                    % In the case of a cap we will consider the volume
                    % element directly in trigauss_cap.
                    [A,B]=meshgrid(yy(1,1:n+1,1),yy(2,1:((n+1)/2),1));
                    yys(:,1)=A(:);
                    yys(:,2)=B(:);
                    [C,D]=meshgrid(yy(1,1:n+1,2),yy(2,1:((n+1)/2),2));
                    yys(:,3)=C(:).*D(:);
                else
                    yy=ones(d,n+2,2);
                    % We need even number of nodes in [0,2*pi] to have a
                    % central simmetry.
                    % yy(1,1:n+2,:)=trigauss_2012(n+1,x(1,1),x(1,2));
                    yy(1,1:n+2,:)=circle_trigquad(n+1,x(1,1),x(1,2));
                    yy(2,1:n+1,:)=trigauss_abssin(n,x(2,2));

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

            else % SPHERICAL QUADRANGLE.
                % fprintf('\n \t SPHERICAL QUADRANGLE');
                yy=ones(d,n+2,2);
                yy(1,1:n+1,:)=trigauss_2012(n,x(1,1),x(1,2));
                yy(2,1:n+2,:)=trigauss_2012(n+1,x(2,1),x(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

            % fprintf('\n \t SPHERICAL QUADRANGLE CC');
            [x1,w1]=trig_cc(2*n+1,x(1,1),x(1,2),method);
            [x2,w2]=trig_cc(2*n+2,x(2,1),x(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


% Convert all points in cartesian coordinates:
% phi=thetaw(:,1), theta=thetaw(:,2), w=thetaw(:,3)
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.

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


if fig==1
    scatter_plot = xw(:,1:3);

    [u,v]=meshgrid(x(1,1):(x(1,2)-x(1,1))/20:x(1,2),x(2,1):(x(2,2)-x(2,1))/20:x(2,2));
    x=(sin(v)).*cos(u);
    y=(sin(v)).*sin(u);
    z=cos(v);
    hold on;
    colormap(0.5*ones(size(z,1),3));
    axis("off");
    mesh(x,y,z)
    plot3(scatter_plot(:,1),scatter_plot(:,2),scatter_plot(:,3),'r.');
    view(180,90);
    hold off;

end



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