
function [wam_cart,wam_sph]=wam_sphrect(x,n)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% We construct a WAM for a spherical rectangle of longitude-latitude type
% of the unit-sphere S^2.
%
% The mesh will be described
%
% 1. in cartesian coordinates ia the N x 3 matrix  "wam_cart";
%
% 2. in spherical coordinates
%
%              x_1=cos(theta_2)*sin(theta_1)
%              x_2=sin(theta_2)*sin(theta_1)
%              x_3=cos(theta_1)
%
%    with "theta_1" in "[0,pi]" and "theta_2" in "[0,2*pi]", via the N x 2
%    matrix "wam_sph".
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% x:  vector 2 x 2, determining the spherical rectangle over which the WAM
%     is computed.
%
%     In particular
%         x(1,:)=[theta1(1) theta1(2)]   x(2,:)=[theta2(1) theta2(2)]
%
%     Note that
%          [theta1(1),theta1(2)] is a subinterval of [0,pi],
%          [theta2(1),theta2(2)] is a subinterval of [0,2*pi].
%
% n: degree of the WAM.
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% wam_cart: it is a (2n+1)^2 x 3 matrix describing the WAM in cartesian
%          coordinates;
%
% wam_sph: it is a (2n+1)^2 x 2 matrix describing the WAM in spherical
%          coordinates.
%--------------------------------------------------------------------------
% Routines:
%--------------------------------------------------------------------------
% All the needed routines are contained in this file (see attached files).
%--------------------------------------------------------------------------
% Note:
%--------------------------------------------------------------------------
% a) The algorithm is optimized for caps.
%
% b) If the region is a cap with the centre in the North Pole, then the WAM
% will have at most (2n+1)*(n+1) - n points.
%--------------------------------------------------------------------------
% Example.
%--------------------------------------------------------------------------
% >> % DEGREE 1 WAM ON THE CAP [0 pi/3; 0 2*pi].
% >> [wam_cart,wam_sph]=wam_sphrect([0 pi/3; 0 2*pi],1)
% 
% wam_cart =
% 
%    -0.3903   -0.6760    0.6250
%    -0.3903    0.6760    0.6250
%    -0.0000   -0.0000    1.0000
%     0.7806         0    0.6250
% 
% 
% wam_sph =
% 
%     0.8957    4.1888
%     0.8957    2.0944
%     0.0000    4.1888
%     0.8957         0
% 
% [wam_cart,wam_sph]=wam_sphrect([0 pi/3; 0 pi/4],1)
% 
% wam_cart =
% 
%     0.6159    0.5536    0.5606
%     0.7651    0.3169    0.5606
%     0.8269    0.0440    0.5606
%     0.3718    0.3343    0.8660
%     0.4619    0.1913    0.8660
%     0.4993    0.0266    0.8660
%     0.0531    0.0477    0.9974
%     0.0660    0.0273    0.9974
%     0.0713    0.0038    0.9974
% 
% 
% wam_sph =
% 
%     0.9757    0.7322
%     0.9757    0.3927
%     0.9757    0.0532
%     0.5236    0.7322
%     0.5236    0.3927
%     0.5236    0.0532
%     0.0715    0.7322
%     0.0715    0.3927
%     0.0715    0.0532
% 

% >>
%--------------------------------------------------------------------------
%%
%% 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@math.unipd.it>
%%          Marco Vianello   <marcov@math.unipd.it>
%%
%% Date    : May 15, 2013
%% Modified: January 4, 2026
%--------------------------------------------------------------------------

if ((x(2,2)-x(2,1)==pi) && x(1,1)==-x(1,2))

    [wam_sph,wam_cart]=wam_cap(x,n); % CAP.

else
    if (x(2,2)-x(2,1)==2*pi) && x(1,1)==0 % CAP.

        x=[-x(1,2) x(1,2); 0 pi];
        [wam_sph,wam_cart]=wam_cap(x,n);

    else % SPHERICAL RECTANGLE.

        [wam_sph,wam_cart]=wam_sph_rect(x,n);

    end
end









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

%--------------------------------------------------------------------------
% wam_cap
%--------------------------------------------------------------------------

function [wam_sph,wam_cart]=wam_cap(x,n)

%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% x = vector 2 x 2, determining the spherical cap over which the WAM
%     is computed.
%     In each row, x must denote the start point and the end point.
%     Remember that the first row elements must be contained in [0,pi],
%     while the second row elements must be in [0,2*pi]).
%     For example, [a b; c d] is a rectangular domain with the first
%     variable theta_1 in [a,b] and the second variable theta_2 in [c,d].
%
% n = degree of the WAM.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% wam_sph: it is a O(2n^2) x 2 matrix describing the WAM in
%          spherical coordinates.
% 
% wam_cart: it is a O(2n^2) x 3 matrix describing the WAM in
%          cartesian coordinates;
%--------------------------------------------------------------------------
% NEEDED ROUTINES.
%--------------------------------------------------------------------------
% 1. wam_arc
% 2. wam_circle
% 3. sph_2_cart
% They are all attached to this file.
%--------------------------------------------------------------------------

wam_polar1=wam_arc(x(1,:),n);
wam_polar1=wam_polar1(1:n+1);
wam_polar2=wam_circle(n);
[wam_polar_th1,wam_polar_th2]=meshgrid(wam_polar1,wam_polar2);
wam_sph=[wam_polar_th1(:) wam_polar_th2(:)];

if nargout > 1
    wam_cart=sph_2_cart(wam_sph);
    [wam_cart,irows]=uniquetol(wam_cart,10^(-15),'ByRows',true);
    wam_sph=wam_sph(irows,:);
end



%--------------------------------------------------------------------------
% wam_sph_rect
%--------------------------------------------------------------------------

function [wam_sph,wam_cart]=wam_sph_rect(x,n)

%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% x : vector 2 x 2, determining the spherical rectangle over which the WAM
%     is computed.
%     In each row, x must denote the start point and the end point.
%     Remember that the first row elements must be contained in [0,pi],
%     while the second row elements must be in [0,2*pi]).
%     For example, [a b; c d] is a rectangular domain with the first
%     variable theta_1 in [a,b] and the second variable theta_2 in [c,d].
%
% n : degree of the WAM.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% wam_sph  : (2n+1)^2 x 2 matrix describing WAM in spherical coordinates;
% wam_cart : (2n+1)^2 x 3 matrix describing WAM in cartesian coordinates.
%--------------------------------------------------------------------------
% NEEDED ROUTINES.
%--------------------------------------------------------------------------
% 1. wam_arc
% 2. sph_2_cart
% The routines are attached to this file.
%--------------------------------------------------------------------------

wam_polar1=wam_arc(x(1,:),n);
wam_polar2=wam_arc(x(2,:),n);
[wam_polar_th1,wam_polar_th2]=meshgrid(wam_polar1,wam_polar2);
wam_sph=[wam_polar_th1(:) wam_polar_th2(:)];

if nargout > 1
    wam_cart=sph_2_cart(wam_sph);
end




%--------------------------------------------------------------------------
% wam_arc
%--------------------------------------------------------------------------

function [wam_polar,wam_cart]=wam_arc(x,n)

%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% x = vector 2 x 1, determining the interval of the unit circle over which
%     the WAM is computed, i.e. x denotes its start point and the end
%     point.
%
% n = degree of the WAM.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% wam_polar: WAM in polar coordinates.
% wam_cart : WAM in cartesian coordinates.
%--------------------------------------------------------------------------
% NEEDED ROUTINES.
%--------------------------------------------------------------------------
% 1. wam_arc
% 2. sph_2_cart
% The routines are attached to this file.
%--------------------------------------------------------------------------

aa=x(1);
bb=x(2);
omega=(bb-aa)/2;
alpha=sin((omega)/2);

j=1:2*n+1; j=j';
tau=cos((2*j-1)*pi/(2*(2*n+1)));
wam_polar=aa+omega+2*asin(alpha*tau);

if nargout > 1
    wam_cart=[cos(wam_polar) sin(wam_polar)];
end




%--------------------------------------------------------------------------
% wam_circle
%--------------------------------------------------------------------------

function [wam_polar,wam_cart]=wam_circle(n)

%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% n = degree of the WAM of the unit circle.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% wam_polar: WAM in polar coordinates.
% wam_cart : WAM in cartesian coordinates.
%--------------------------------------------------------------------------
% NEEDED ROUTINES.
%--------------------------------------------------------------------------
% 1. . sph_2_cart
% The routine is attached to this file.
%--------------------------------------------------------------------------

wam_polar=linspace(0,2*pi,2*n+2);
wam_polar=(wam_polar(1:end-1))';

if nargout > 1
    wam_cart=[cos(wam_polar) sin(wam_polar)];
end





%--------------------------------------------------------------------------
% sph_2_cart
%--------------------------------------------------------------------------

function pts=sph_2_cart(theta)

%--------------------------------------------------------------------------
% INPUT.
%--------------------------------------------------------------------------
% theta = polar coordinates of the set.
%--------------------------------------------------------------------------
% OUTPUT.
%--------------------------------------------------------------------------
% pts = conversion from polar to cartesian coordinates of the set.
%--------------------------------------------------------------------------
% NEEDED ROUTINES.
%--------------------------------------------------------------------------
% No routine is needed.
%--------------------------------------------------------------------------

theta_1=theta(:,1);
theta_2=theta(:,2);

pts(:,1)=cos(theta_2).*sin(theta_1);
pts(:,2)=sin(theta_2).*sin(theta_1);
pts(:,3)=cos(theta_1);
