% Random test for cubature rule
% Testing polynomial exactness for small degrees, 

format long e;

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

exactness_deg = 3;

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

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

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

% Sympy result
for i=1:size(monomials,1)
    poli = monomials(i,:);
    command = sprintf( ...
'python integrate-python.py "x**%d*y**%d*z**%d" %0.50f %0.50f %0.50f %0.50f' ...
    ,poli(1),poli(2),poli(3),region(1,1),region(1,2),region(2,1),region(2,2));
    [ret,sympy_res] = system(command);
    sympy_res = str2num(sympy_res);

%    disp(sympy_res)

    % Cubature formula result
    [tw,xw] = cubature_2_torus(region,exactness_deg,0);

    %Cubature result
    x = xw(:,1);
    y = xw(:,2);
    z = xw(:,3);
    w = xw(:,4);
    summands = x.^(poli(1)).*y.^(poli(2)).*z.^(poli(3)).*w;

    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' ...
    ,poli(1),poli(2),poli(3),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);

