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

% this program is optimized for caps in S^2
% for general S^d use cubature_d_sphere
%
% let use the (iper)spherical coordinate system over S^2 (in R^3)
%
% x_1=cos(phi)sin(theta)
% x_2=sin(phi)sin(theta)
% x_3=cos(theta)
%
% with phi in [0,2*pi] and theta in [0,pi]
%
% remember the volume element for these coordinates
% sin(theta) (formally abs(sin(theta)))
%
% input
% x = vector 2 x 2
% horizontally, x must denote the start point and the end point for each
% range [phi_s,phi_f], [theta_s,theta_f] (remember that the first
% one must be contained in [0,2*pi] and the other in [0,pi])
%
% n = degree of polinomially exactness of nodes and weights in output
%
% fig = 0,1
% if 1
% we will use a bidimensional Delaunay triangulation over the first two
% coordinates of xw the color will be defined with the average of the
% weights of the three vertices of each triangle
%
% output
% thetaw = vector (n+2)*(n+1) x 3
% for each row it will denote the iperspherical coordinates of nodes
% and in the last column it will denote the weights for each node
%
% xw = vector (n+2)*(n+1) x 4
% for each row it will denote the cartesian coordinates of nodes
% in the last column it will denote the weights for each node
%
% if x represents a cap, we will have less number of nodes, that is
% if n is odd
% thetaw = vector (n+1)*(n+1)/2 x 3
% xw     = vector (n+1)*(n+1)/2 x 4
% if n is even
% thetaw = vector (n+1)*(n+2)/2-(n/2) x 3
% xw     = vector (n+1)*(n+2)/2-(n/2) x 4
%
% functions not usually lie in standard Matlab
%
% trigauss.m
% from Prof.Vianello web page from Padua University
% http://www.math.unipd.it/~marcov/CAAsoft.html
%
% trigauss_sph.m
%
% pick.m
% from the web site www.mathworks.com
% http://www.mathworks.com/matlabcentral/fileexchange/12724-picking-element
% s-from-a-set/content/pick.m
tic
d=size(x,1);
if d~=2
    fprintf('Use cubature_d_sphere\n')
else

if (x(1,2)-x(1,1)==2*pi) && (x(2,1)==0)
    if mod(n,2)==1
        yy=ones(d,n+1,2);
        yy(1,1:n+1,:)=trigauss(n,x(1,1),x(1,2));
        yy(2,1:n+1,:)=trigauss_cap(n,x(2,2));
        % note that in the case of a cap we will consider the volume
        % element directly in trigauss_sph, so we will not need to add it
        [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(:);
    end
    if mod(n,2)==0
        yy=ones(d,n+2,2);
        % we will need even number of nodes in [0,2*pi] to have a
        % a central simmetry, so we will use the degree n+1
        yy(1,1:n+2,:)=trigauss(n+1,x(1,1),x(1,2));
        yy(2,1:n+1,:)=trigauss_cap(n,x(2,2));
        % note that in the case of a cap we will consider the volume
        % element directly in trigauss_sph, so we will not need to add it
        [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

yy=ones(d,n+2,2);

fprintf('ff\n')

yy(1,1:n+1,:)=trigauss(n,x(1,1),x(1,2));
yy(2,1:n+2,:)=trigauss(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
% convert all points in cartesian coordinates
xw=ones(size(thetaw,1),d+2);
xw(:,d+2)=thetaw(:,d+1);
xw(:,1)=cos(thetaw(:,1));
xw(:,2)=sin(thetaw(:,1));
xw(:,1)=xw(:,1).*sin(thetaw(:,2));
xw(:,2)=xw(:,2).*sin(thetaw(:,2));
xw(:,d+1)=cos(thetaw(:,d));


% figures

if size(x,1)==2
if fig==1
        pos=find(xw(:,3)>=0);
        neg=find(xw(:,3)<0);
        destra=find(xw(:,2)>=0);
        sinistra=find(xw(:,2)<0);
        anteriore=find(xw(:,1)>=0);
        posteriore=find(xw(:,1)<0);
        if size(pos,1)~=0
           xwpos=xw(pos,:);
           tri1=delaunay(xwpos(:,1),xwpos(:,2));
           colore1=(xwpos(tri1(:,1),4)+xwpos(tri1(:,2),4)+xwpos(tri1(:,3),4))/3;
        end
        if size(neg,1)~=0
           xwneg=xw(neg,:);
           tri2=delaunay(xwneg(:,1),xwneg(:,2));
           colore2=(xwneg(tri2(:,1),4)+xwneg(tri2(:,2),4)+xwneg(tri2(:,3),4))/3;
        end
        if size(destra,1)~=0
           xwd=xw(destra,:);
           trid=delaunay(xwd(:,1),xwd(:,3));
           colored=(xwd(trid(:,1),4)+xwd(trid(:,2),4)+xwd(trid(:,3),4))/3;
        end
        if size(sinistra,1)~=0
           xws=xw(sinistra,:);
           tris=delaunay(xws(:,1),xws(:,3));
           colores=(xws(tris(:,1),4)+xws(tris(:,2),4)+xws(tris(:,3),4))/3;
        end
        if size(anteriore,1)~=0
           xwa=xw(anteriore,:);
           tria=delaunay(xwa(:,2),xwa(:,3));
           colorea=(xwa(tria(:,1),4)+xwa(tria(:,2),4)+xwa(tria(:,3),4))/3;
        end
        if size(posteriore,1)~=0
           xwp=xw(posteriore,:);
           trip=delaunay(xwp(:,2),xwp(:,3));
           colorep=(xwp(trip(:,1),4)+xwp(trip(:,2),4)+xwp(trip(:,3),4))/3;
        end
        figure(1)
        hold on
        if size(pos,1)~=0
           trisurf(tri1,xwpos(:,1),xwpos(:,2),xwpos(:,3),colore1)
        end
        if size(neg,1)~=0
           trisurf(tri2,xwneg(:,1),xwneg(:,2),xwneg(:,3),colore2)
        end
        if size(destra,1)~=0
           trisurf(trid,xwd(:,1),xwd(:,2),xwd(:,3),colored)
        end
        if size(sinistra,1)~=0
           trisurf(tris,xws(:,1),xws(:,2),xws(:,3),colores)
        end
        if size(anteriore,1)~=0
           trisurf(tria,xwa(:,1),xwa(:,2),xwa(:,3),colorea)
        end
        if size(posteriore,1)~=0
           trisurf(trip,xwp(:,1),xwp(:,2),xwp(:,3),colorep)
        end
        hold off
    end
end

end
toc