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

% this program is NOT optimized for caps in S^2
%
% let use the (hyper)spherical coordinate system over S^d (in R^(d+1))
%
% x_1=cos(phi)sin(theta_1)....sin(theta_(d-1))
% x_2=sin(phi)sin(theta_1)....sin(theta_(d-1))
% x_3=cos(theta_1)sin(theta_2)....sin(theta_(d-1))
% x_4=cos(theta_2)sin(theta_3)....sin(theta_(d-1))
% ....
% x_d=cos(theta_(d-2))sin(theta_(d-1))
% x_(d+1)=cos(theta_(d-1))
%
% with phi in [0,2*pi] and theta_i in [0,pi]
%
% remember the volume element for these coordinates
% sen(theta_(d-1))^(d-1)sen(theta_(d-2))^(d-2)....sen(theta_1)
%
% input
% x = vector d x 2
% first dimension of x will denote the dimension of the space R^(d+1)
% where we are working: if size(x,1)=d then we will consider S^d;
% horizontally, x must denote the start point and the end point for each
% range [phi_s,phi_f], [theta_i_s,theta_i_f] (remember that the first
% one must be contained in [0,2*pi] and the others in [0,pi])
%
% n = degree of polinomially exactness of nodes and weights in output
%
% fig = 0,1
% if 1
% if size(x,1)=1 draw the nodes using plot
% if size(x,1)=2 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
% if size(x,1)>=3 do nothing
%
% output
% thetaw = vector (n+d)...(n+1) x (d+1)
% 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+d)...(n+1) x (d+2)
% for each row it will denote the cartesian coordinates of nodes
% in the last column it will denote the weights for each node
%
% 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
%
% 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);

yy=ones(d,n+1+d-1,2);
yy=-2.*yy;
% initialize the matrix with -2 elements: it will be useful later

for i=1:1:d % let use trigauss code for each dimension
    yy(i,1:n+1+(i-1),:)=trigauss(n+(i-1),x(i,1),x(i,2));
end
% in the degree of trigauss we will consider also the volume element

% yy(i,:,1) contain the nodes of each interval
% yy(i,:,2) contains the weights of each interval

xx=pick(1:n+d,d,'or');
% alternative method for ndgrid: it will obtain all the combinations
% of index of yy so we will have all the coordinates of points in
% the tensor product space

% use xx as index in yy
yys=ones(size(xx,1),d+1);
for j=1:1:d
    yys(:,j)=yy(j,xx(:,j),1);
    yys(:,d+1)=yys(:,d+1).*yy(j,xx(:,j),2)';
end


% we have gained all the coordinates of all points more some
% unnecessary points (in trigauss we used a matrix with n+d columns
% but we needed only n+i of them)
% let eliminate this unnecessary points (it will be useful to have
% initialized the yy with -2 elements)
cont=0;
mm=min(yys(:,1:d),[],2);
for i=1:1:length(mm)
    if mm(i)>-1
        cont=cont+1;
        ind_buono(cont)=i;
    end
end

thetaw=yys(ind_buono,:);

% scale the weights using the jacobian matrix of iperspherical coordinates
% sen(theta_(d-1))^(d-1)sen(theta_(d-2))^(d-2)....sen(theta_1)
if d>=2 
    for j=d:-1:2
        thetaw(:,d+1)=thetaw(:,d+1).*((sin(thetaw(:,j)).^(j-1)));
    end
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));
if d>=2
   for i=2:1:d
       xw(:,1)=xw(:,1).*sin(thetaw(:,i));
       xw(:,2)=xw(:,2).*sin(thetaw(:,i));
   end
   for j=3:1:d
       xw(:,j)=cos(thetaw(:,j-1));
       for i=j+1:1:d+1
           xw(:,j)=xw(:,j).*sin(thetaw(:,i-1));
       end
   end
   xw(:,d+1)=cos(thetaw(:,d));
end

% figures
if size(x,1)==1
    if fig==1
        figure(1)
        plot(xw(:,1),xw(:,2),'b*');
    end
end

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
toc