% Random test for cubature rule
% Testing polynomial exactness for small degrees, 
degs = [2,3,4,5,10];

output_hist = [];
for jj=1:length(degs)

format long e;

% number of tests
abs_err_hist = [];
rel_err_hist = [];
val_hist = [];

exactness_deg = degs(jj);

monomials = [];

for i = 0:exactness_deg
    for j = 0:exactness_deg-i
        for k=0:exactness_deg-(i+j)
        monomials = [monomials; i,j,k];
        end
    end
end
disp('Numero di integrali:');
size(monomials,1)

region = [pi/6,pi/3;0,pi/2];

% "Standard" torus with radii 2,1
R = [2,1];
flag = 0;

[tw,xw] = cubature_2_torus(region,exactness_deg);
% Sympy result
for i=1:size(monomials,1)
    px = monomials(i,1);
    py = monomials(i,2);
    pz = monomials(i,3);

    polit = @(t,p) (2+cos(t)).^(px+py+1).*(cos(p)).^(px).*(sin(p)).^(py).*(sin(t)).^(pz);
    sympy_res = dblquad(polit,region(1,1),region(1,2),region(2,1),region(2,2));

    %Cubature result
    x = xw(:,1);
    y = xw(:,2);
    z = xw(:,3);
    w = xw(:,4);
    summands = x.^(px).*y.^(py).*z.^(pz).*w;
    % Unused polar calculations
    %fx = feval(@(t,p) (2+cos(t)).^(px+py).*(cos(p)).^(px).*(sin(p)).^(py).*(sin(t)).^(pz),tw(:,1),tw(:,2));
    %summands = fx.*tw(:,3);

    quad_result = sum(summands);
    val_hist = [val_hist; quad_result];
    if ((sympy_res == 0) || (quad_result == 0))
        abs_err_hist = [abs_err_hist; abs(sympy_res - quad_result)];
        rel_err_hist = [rel_err_hist; abs(sympy_res - quad_result)];
    else
        abs_err_hist = [abs_err_hist; abs(sympy_res - quad_result)];
        rel_err_hist = [rel_err_hist; ...
            abs(abs(sympy_res - quad_result)/quad_result)];
    end
    
    output = sprintf( ...
'Integrato il poliomio x^%dy^%dz^%d con risultato %1.30e ed errore %1.30e' ...
    ,px,py,pz,quad_result,abs(sympy_res - quad_result));
    disp(output);

%    if ((abs(sympy_res - sum(summands))) > 10^(-1))
%		disp('Test1: Approximation not good at degree')
%		ni
%        flag = 1;
%    end
end

bool = flag

final_output = sprintf( ...
'La formula di grado %d integra esattamente tutti i polinomi di grado <= %d a meno di:\n\tErrore assoluto massimo: %1.8e\n\tErrore relativo massimo: %1.8e' ...
, exactness_deg, exactness_deg, max(abs_err_hist),max(rel_err_hist));
disp(final_output);

output_hist = [output_hist; exactness_deg, size(monomials,1), max(abs_err_hist),max(rel_err_hist)];

end

for i=1:size(output_hist,1)
    printf('%d & %d & %1.5e & %1.5e \\\\ ',output_hist(i,1),output_hist(i,2),output_hist(i,3),output_hist(i,4));
end
