function [th,y]=wam_s2(x,deg,fig)
% We will cosntruct a WAM for a subregion of a 2-dimensional sphere S^2
% let use the 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]
%
% 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 others in [0,pi])
%
% deg = degree of the WAM
%
% fig = 0,1
% fig=1 create a figure using plot3
%
% output
% th = (2n+1)^2 x 2
% WAM in spherical coordinates
%
% y = (2n+1)^2 x 3
% WAM in cartesian coordinates
%
%
% the algorithm is optimized for caps:
% if the region is a cap with the centre in the North Pole, then the WAM
% will have half points

if (x(1,2)-x(1,1)==2*pi) && x(2,1)==0
    alpha=sin(x(2,2)/2);
    theta2=2.*asin(alpha.*cos(((2.*(1:1:deg)-1).*pi)./(2.*(2.*deg+1))));
    theta1=2.*pi.*(0:1:2*deg+1)./(2*deg+2);

else

    for i=1:1:length(x(:,1))
        xx1(i,1)=0;
        xx1(i,2)=(x(i,2)-x(i,1));
    end
    for i=1:1:length(x(:,1))
        %xx2(i,1)=-(xx1(i,1)+xx1(i,2))/2;
        xx2(i,2)=(xx1(i,1)+xx1(i,2))/2;
        alpha(i)=sin(xx2(i,2)/2);
    end

    theta1=2.*asin(alpha(1).*cos(((2.*(1:1:2*deg+1)-1).*pi)./(2.*(2.*deg+1))));
    theta2=2.*asin(alpha(2).*cos(((2.*(1:1:2*deg+1)-1).*pi)./(2.*(2.*deg+1))));

    theta2(:)=theta2(:)+(xx1(2,1)+xx1(2,2))/2;
    theta2(:)=theta2(:)+x(2,1);
    theta1(:)=theta1(:)+(xx1(1,1)+xx1(1,2))/2;
    theta1(:)=theta1(:)+x(1,1);


end

[A,B]=meshgrid(theta1(:),theta2(:));
th(:,1)=A(:);
th(:,2)=B(:);
if (x(1,2)-x(1,1)==2*pi) && x(2,1)==0
    th(end+1,2)=0;
end

y(:,1)=sin(th(:,2)).*cos(th(:,1));
y(:,2)=sin(th(:,2)).*sin(th(:,1));
y(:,3)=cos(th(:,2));

if fig
    figure(1)
    plot3(y(:,1),y(:,2),y(:,3),'o')
end