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

format long e;

%%%%%
%
% INPUT:
%   Required arguments:
%   x:      Geographical rectangle in toroidal coordinates
%   n:      Degree of precision
%    
%   Optional arguments:
%   fig:    Boolean; plot points if true (>0)
%   method: 
%       1 for Gaussian points
%       2 for Antigaussian points
%       3 for Gauss-Kronrod points
%       (note: if method=2,3, then
%       "n" refers to the degree of 
%       the basic gaussian formula)
%   radii:  Torus radii
%
% OUTPUT:
%   thetaw: nodes (in toroidal coordinates) and weights
%   xw: nodes (in cartesian coordinates) and weights
%
%%%%%

% Set some defaults for the optional arguments
if nargin < 5
    radii = [2,1];
end
if nargin < 4
    method = 1;
end
if nargin < 3
    fig = 0;
end

% Check if the input domain has the right dimensions
d=size(x,1);
if d~=2
    fprintf('Wrong size for the domain\n')
else

% Torus radii
R = radii(1);
r = radii(2);

% Init yy to have the right (maximum) dimensions
yy=ones(d,2*n+3,2);

% Calculating angles and weights 
% via trigonometric gaussian formulas
% using the method specified in input argument
vert_extra = 1;
if method == 1
    yy(1,1:n+2,:)=trigauss(n+1,x(1,1),x(1,2));
    yy(2,1:n+1,:)=trigauss(n,x(2,1),x(2,2));
elseif method == 2
    yy(1,1:n+3,:)=antigauss(n+1,x(1,1),x(1,2));
    yy(2,1:n+2,:)=antigauss(n,x(2,1),x(2,2));
    n = n+1;
elseif method == 3
    vert_extra = 2;
    yy(1,1:2*n+3,:)=trigauss_kronrod(n+1,x(1,1),x(1,2));
    yy(2,1:2*n+1,:)=trigauss_kronrod(n,x(2,1),x(2,2));
    n = 2*n;
end

[A,B]=meshgrid(yy(1,1:n+1+vert_extra,1),yy(2,1:n+1,1));

yys(:,1)=A(:);
yys(:,2)=B(:);

[C,D]=meshgrid(yy(1,1:n+1+vert_extra,2),yy(2,1:n+1,2));
% Note that angles are stored in yys(:,1:2) 
% and weights in yys(:,3) 
element = r*(R+ r*cos(yys(:,1)));
yys(:,3)=C(:).*D(:).*element;

thetaw=yys;

% convert all points in cartesian coordinates
% Toroidal coords:
% x_1=(R + r cos(theta))cos(phi)
% x_2=(R + r cos(theta))sin(phi)
% x_3=sin(theta)

xw=ones(size(thetaw,1),d+2);
% Copy the weights as-is
xw(:,d+2)=thetaw(:,d+1);

% xw(:,1:2) are x and y coordinates
xw(:,1)=cos(thetaw(:,2));
xw(:,2)=sin(thetaw(:,2));
% Multiplying for (R + rcos(theta))
xw(:,1)=xw(:,1).*(R + r*cos(thetaw(:,1)));
xw(:,2)=xw(:,2).*(R + r*cos(thetaw(:,1)));
% xw(:,3) is the z-th coordinate
xw(:,d+1)=r*sin(thetaw(:,1));

% figures, if required

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

    [v,u]=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=(2+cos(v)).*cos(u);
    y=(2+cos(v)).*sin(u);
    z=sin(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

end
%toc

% Helper functions

function tw=double_circle_trigquad(n,alpha,beta)

if nargin == 1
    alpha=0;
    beta=4*pi;
end
N=n+1;
w=(4*pi/N)*ones(N,1);
t=linspace(2*pi/N,4*pi-2*pi/N,N); t=alpha+t';

tw=[t w];
